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ABSTRACT 


This  annual  report  describes  our  progress  during  the  period  from  May  1996  to  April  1997. 
Two  tasks  are  described  in  this  report.  The  first  task  considers  the  hybridization  of  the 
finite-element  method  (FEM)  and  the  shooting-and-bouncing-ray  (SBR)  method  for 
scattering  by  large  bodies  with  small,  inhomogeneous  protruding  scatterers,  and  the 
hybridization  of  the  method  of  moments  (MoM)  and  SBR  method  for  scattering  from 
conformal  slotted  waveguide  arrays  on  a  large,  complex  platform.  The  second  task  studies 
a  variety  of  finite-element  and  boundary-integral  (FE-BI)  methods  for  three-dimensional 
electromagnetic  analysis.  Several  journal  articles  and  conference  papers  supported  by  the 
research  are  also  listed  in  this  report. 
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PROGRESS 


The  main  goal  of  our  proposed  research  is  to  develop  numerical  methods  that  can  simulate 
the  interaction  of  electromagnetic  fields  with  complex  systems  in  a  large,  complex 
environment.  During  the  second  year  of  the  award,  we  have  made  significant  progress 
toward  this  goal.  We  have  worked  on  several  tasks,  which  form  the  integral  parts  of  the 
proposed  research  and  whose  successful  completion  is  vital  to  the  development  of  the 
proposed  technique.  Two  of  these  tasks  are  described  in  detail  in  this  report. 

The  first  task  deals  with  the  development  of  two  hybrid  techniques  to  solve  two  different 
classes  of  electromagnetic  scattering  problems.  The  first  technique  combines  the  finite- 
element  method  (FEM)  and  the  shooting-and-bouncing-ray  (SBR)  method  for  scattering  by 
large  bodies  with  small,  inhomogeneous  protruding  scatterers.  The  second  technique 
combines  the  method  of  moments  (MoM)  and  SBR  method  to  characterize  the  scattering  of 
conformal  slotted  waveguide  antenna  arrays  in  a  large,  complex  platform.  This  task  can  be 
considered  as  the  further  extension  of  the  work  described  in  Tasks  1  and  2  in  the  first 
year's  annual  report.  It  demonstrates  further  the  power  of  the  proposed  hybrid  technique 
for  complex  electromagnetic  problems. 

The  second  task  studies  a  variety  of  finite-element  and  boundary-integral  (FE-BI)  methods 
for  three-dimensional  electromagnetic  analysis.  A  new  formulation  is  proposed,  which  is 
shown  to  be  accurate,  efficient,  and  immune  to  interior  resonance  corruption.  This  work 
lays  a  foundation  for  the  development  of  a  very  useful  and  powerful  technique,  which 
combines  FEM  with  the  multilevel  fast  multipole  method  (MLFMM),  for  a  large  class  of 
electromagnetic  problems. 

In  addition  to  the  two  desks  described  above,  we  have  also  (i)  developed  the  concept  of 
complementary  perfectly  matched  layers  (PML)  to  significantly  improve  the  accuracy  of  the 
FEM  solution  (Appendix  5);  (ii)  developed  the  spectral  Lanczos  decomposition  method 
(SLDM)  for  efficient  time-domain  and  frequency-domain  FEM  solution  of  Maxwell's 
equations;  (iii)  applied  the  method  of  asymptotic  waveform  evaluation  (AWE)  for  multi¬ 
frequency  scattering  analysis;  and  (iv)  developed  a  new  k-space  method  using  the 
transpose-free  quasi  minimum  residual  (TFQMR)  method  and  FFT  for  solving  volume- 
integral  equation  arising  from  the  problem  of  electromagnetic  interaction  with 
inhomogeneous  objects  such  as  the  human  body  (Appendix  6). 
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Task  1:  Hybrid  SBR/FEM  and  SBR/MoM  methods  for  large  and  complex 
scattering  problems 

Two  years  ago,  our  research  group  developed  a  hybrid  technique  for  computing  scattering 
by  large  bodies  with  cracks  and  cavities.  This  technique  employs  the  shooting-and- 
bouncing-ray  (SBR)  method  to  compute  the  scattering  by  large  bodies  and  uses  the  finite- 
element  method  (FEM)  to  characterize  the  cracks  and  cavities.  The  two  methods  are 
combined  through  a  coupling  scheme  based  on  the  electromagnetic  equivalence  principle 
and  the  reciprocity  theorem.  The  coupling  scheme  is  designed  in  such  a  manner  that  it 
includes  all  significant  interactions  between  the  FEM  and  the  SBR  method  and  it  permits 
the  SBR  and  FEM  computations  to  be  done  separately.  The  resulting  technique  is  shown  to 
be  efficient  and  accurate.  During  the  first  year  of  the  award,  we  extended  this  method  to  the 
calculation  of  the  radiation  patterns  of  conformal  antennas  in  a  complex  environment  (see 
Task  1  of  the  first  year's  annual  report).  Comparison  with  experimental  data  showed  that 
the  technique  can  predict  accurately  the  effect  of  the  environment  on  the  radiation  of 
conformal  antennas. 

During  the  past  year,  we  further  developed  a  hybrid  FEM/SBR  method  to  compute 
scattering  by  large  bodies  with  small  inhomogeneous  protruding  scatterers,  a  problem  that 
is  very  important,  but  fundamentally  different  from  those  treated  before.  To  be  more 
specific,  we  first  employ  the  field  equivalence  principle  to  replace  the  protruding  scatterers 
by  a  set  of  equivalent  electric  and  magnetic  currents.  The  total  scattered  field  then  becomes 
the  superposition  of  the  field  scattered  by  the  large  body  without  protrusions,  which  is 
calculated  using  the  SBR  method,  and  the  field  radiated  by  the  equivalent  currents  in  the 
presence  of  the  large  body,  which  is  also  calculated  using  the  SBR  method  with  the  aid  of 
the  reciprocity  theorem.  The  required  equivalent  currents  are  computed  using  the  FEM, 
which  permits  the  handling  of  complex  material  composition  of  the  protrusions.  The 
method  has  been  applied  to  two-dimensional  problems,  simply  to  demonstrate  the 
feasibility  of  the  proposed  technique.  The  technical  details  are  given  in  Appendix  1,  which 
includes  some  very  encouraging  results. 

In  the  first  year's  annual  report,  we  reported  our  preliminary  work  on  the  analysis  of 
cylindrically  conformal  waveguide  slot  antennas  and  arrays.  Since  then  we  have  completed 
the  analysis  and  studied  the  scattering  properties  of  a  variety  of  conformal  waveguide  slot 
antennas  and  arrays  (Appendix  2).  Furthermore,  we  have  combined  the  method  of  analysis 
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with  the  SBR  method  to  compute  scattering  from  conformal  waveguide  slot  antennas  and 
arrays  in  a  complex  geometry.  This  extension  is  important  because  almost  all  conformal 
antennas  and  arrays  are  placed  on  a  complex  geometry  in  practice  and  their  scattering 
characterization  is  critical  for  applications  such  as  target  identification.  The  technical  details 
are  given  in  Appendix  3. 


Task  2:  Hybrid  Hnite-element  and  boundary-integral  method  for  electro¬ 
magnetic  analysis 

In  our  originally  proposed  research,  we  planned  to  develop  hybrid  techniques  to  combine  a 
high-frequency  asymptotic  method  with  a  numerical  method  to  simulate  the  interaction  of 
electromagnetic  fields  with  complex  systems  in  a  large,  complex  platform.  To  date,  we 
have  successfully  developed  several  such  techniques  as  described  in  Task  1  of  the  first 
year's  annual  report  and  in  the  Task  1  of  this  report.  Because  of  the  use  of  the  high- 
frequency  method,  the  source  of  the  electromagnetic  fields  must  be  far  away  from  the 
platform,  which  is  the  case  for  many  practical  problems.  However,  there  is  another  class  of 
important  problems  which  involve  sources  near  or  on  the  platform,  such  as  the  problem  of 
mutual  coupling  between  two  antennas  on  the  same  platform  and  the  interaction  of  a 
radiating  source  with  a  system  on  the  same  body.  Such  a  problem  cannot  be  analyzed  by 
the  method  proposed  in  this  research.  We  decided  to  go  beyond  the  originally  proposed 
plan  and  develop  a  technique  for  this  class  of  problems. 

The  technique  which  is  to  be  developed  combines  FEM  with  the  multilevel  fast  multipole 
method  (MLFMM).  It  is  based  on  the  finite-element  and  boundary-integral  (FE-BI)  method 
pioneered  by  the  PI.  The  FE-BI  method  first  divides  the  problem  into  an  interior  and 
exterior  problems.  The  field  in  the  interior  region  is  formulated  using  FEM,  and  the  field  in 
the  exterior  region  is  represented  by  a  boundary-integral  equation  (BIE).  The  interior  and 
exterior  fields  are  then  coupled  by  the  field  continuity  conditions.  Although  the  FE-BI 
method  is  remarkably  more  powerful  than  other  numerical  techniques  in  dealing  with 
complex  problems,  it  still  has  a  bottleneck  which  is  the  full  matrix  generated  by  BIE.  This 
bottleneck  severely  limits  the  capability  of  the  FE-BI  method  in  dealing  with  large 
problems.  Our  objective  is  to  apply  MLFMM  to  BIE  to  completely  remove  the  bottleneck  in 
the  FE-BI  method  for  general  3D  problems.  The  first  step  is  to  study  a  variety  of  FE-BI 
formulations  and  identify  the  most  accurate  and  efficient  one. 
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During  the  past  year,  we  studied  in  detail  a  variety  of  formulations  for  the  hybrid  FE-BI 
method  for  3D  electromagnetic  scattering  by  inhomogeneous  objects.  It  is  shown  that  the 
efficiency  and  accuracy  of  the  FE-BI  method  depend  highly  on  the  formulation  and 
discretization  of  the  BIE  used.  A  simple  analysis  of  matrix  condition  identifies  the 
efficiency  of  the  different  FE-BI  formulations  and  an  analysis  of  weighting  functions 
shows  that  the  traditional  FE-BI  formulations  cannot  produce  accurate  solutions.  A  new 
formulation  is  then  proposed  and  numerical  results  show  that  the  resulting  solution  has  a 
good  efficiency  and  accuracy  and  is  completely  immune  to  the  problem  of  interior 
resonance.  The  technical  details  are  given  in  Appendix  4. 


FUTURE  WORK 


In  the  coming  year,  we  plan  to  focus  our  effort  on  the  development  of  the  FEM/MLFMM 
for  fast  and  accurate  analysis  of  conformal  antennas  and  other  electromagnetic  devices  on  a 
large,  complex  platform.  Such  a  technique  is  urgently  needed  for  many  applications,  such 
as  EMP  due  to  nearby  sources,  design  of  special-purpose  antennas  taking  into 
consideration  of  the  environment/platform,  and  characterization  of  mutual  coupling  between 
two  antennas  on  a  large,  complex  body.  Currently,  MLFMM  has  been  implemented  only 
for  conducting  and  impedance  surfaces,  both  of  which  contain  only  one  type  of  integral 
operator.  Its  application  to  the  proposed  work  requires  the  treatment  of  two  different  types 
of  integral  operators.  In  addition,  we  also  plan  to  develop  fast  solvers  based  on  the  spectral 
Lanczos  decomposition  method  (SLDM)  and  the  method  of  asymptotic  waveform 
evaluation  (AWE)  for  the  FEM  analysis. 
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APPENDIX  1 


Hybrid  FEM/SBR  Method  to  Compute  Scattering 
by  Large  Bodies  with  Small  Protruding  Scatterers 

X.  Q.  Sheng  and  J.  M.  Jin 
Center  for  Computational  Electromagnetics 
Department  of  Electrical  and  Computer  Engineering 
University  of  Illinois  at  Urbana-Champaign 
Urbana,  Illinois  61801-2991, USA 

Abstract-A  hybrid  method  that  combines  the  finite-element  method  (FEM) 
and  the  shooting-and-bouncing-ray  (SBR)  method  is  presented  to  compute  scat¬ 
tering  by  large  bodies  with  small  protruding  scatterers.  In  the  method,  the  field 
equivalence  principle  is  employed  to  replace  the  protruding  scatterers  by  a  set  of 
equivcilent  electric  and  magnetic  currents.  The  total  scattered  field  then  becomes 
the  superposition  of  the  field  scattered  by  the  large  body  without  protrusions, 
which  is  calculated  using  the  SBR  method,  and  the  field  radiated  by  the  equiva¬ 
lent  currents  in  the  presence  of  the  large  body,  which  is  also  calculated  using  the 
SBR  method  with  the  aid  of  the  reciprocity  theorem.  The  required  equivalent 
currents  are  computed  using  the  FEM,  which  permits  the  handling  of  complex 
material  composition  of  the  protrusions.  Two-dimensional  examples  are  presented 
to  demonstrate  the  feasibility  of  the  proposed  method. 
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I.  INTRODUCTION 

Recently,  a  hybrid  technique  was  developed  for  computing  scattering  by  large 
bodies  with  cracks  and  cavities  [1].  This  technique  employs  the  shooting-and- 
bouncing-ray  (SBR)  method  to  compute  the  scattering  by  large  bodies  and  uses 
the  finite-element  method  (FEM)  to  characterize  the  cracks  and  cavities.  The  two 
methods  are  combined  through  a  coupling  scheme  based  on  the  electromagnetic 
equivalence  principle  and  the  reciprocity  theorem.  The  coupling  scheme  is  de¬ 
signed  in  such  a  manner  that  it  includes  all  significant  interactions  between  the 
FEM  and  the  SBR  method  and  it  permits  the  SBR  and  FEM  computations  to 
be  done  separately.  The  resulting  technique  is  shown  to  be  efficient  and  accurate 
and,  because  of  this,  it  is  extended  to  the  calculation  of  the  radiation  patterns  of 
conformal  antennas  in  a  complex  environment  [2]. 

The  hybrid  FEM/SBR  method  summarized  above  belongs  to  a  larger  class  of 
hybrid  method  pioneered  by  Thiele  et  al.  [3]-[ll].  A  common  feature  in  these 
methods  is  that  they  combine  a  high-frequency  asymptotic  technique,  such  as  the 
geometrical  optics  (GO),  physical  optics  (PO),  or  geometrical  theory  of  diffraction 
(GTD),  with  a  low-frequency  numerical  technique,  such  as  the  method  of  moments 
(MoM),  FEM,  or  finite-difference  method  (FDM).  They  are  specifically  designed 
to  tackle  three  types  of  problems  that  cannot  be  handled  accurately  and  efficiently 
by  either  a  high-frequency  method  or  a  low-frequency  method.  The  first  type 
includes  scatterers  with  a  size  in  the  intermediate  region  between  the  high  and 
low  frequencies  [7]-[ll].  The  second  type  is  the  scattering  by  small  objects  in  the 
presence  of  large  bodies  [3]-[6].  The  third  type  is  the  scattering  by  large  bodies 
having  small  indenting  structures  such  as  cracks,  gaps,  and  cavities  [1],  [2]. 

In  this  article,  we  extend  the  hybrid  FEM/SBR  method,  developed  originally 
for  the  third  type  of  problem,  to  compute  scattering  by  large  bodies  with  small 
protruding  scatterers,  a  problem  that  belongs  the  second  type.  To  be  more  specific, 
we  first  employ  the  field  equivalence  principle  to  replace  the  protruding  scatterers 
by  a  set  of  equivalent  electric  and  magnetic  currents.  The  total  scattered  field 
then  becomes  the  superposition  of  the  field  scattered  by  the  large  body  without 
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protrusions,  which  is  calculated  using  the  SBR  method,  and  the  field  radiated  by 
the  equivalent  currents  in  the  presence  of  the  large  body,  which  is  also  calculated 
using  the  SBR  method  with  the  aid  of  the  reciprocity  theorem.  The  required 
equivalent  currents  are  computed  using  the  FEM,  which  permits  the  handling  of 
complex  material  composition  of  the  protrusions.  In  this  article,  the  formulation 
and  analysis  are  described  for  two-dimensional  problems,  simply  to  demonstrate 
the  feasibility  of  the  proposed  technique. 

II.  FORMULATION 

Consider  the  problem  of  wave  scattering  by  a  large,  perfectly  conducting  body 
with  a  small  protruding  structure,  whose  cross-section  is  illustrated  in  Fig.  1.  The 
protruding  structure  can  be  a  perfect  conductor  or  a  dielectric/magnetic  material 
or  a  combination  of  these.  In  accordance  with  the  field  equivalence  principle  [12], 
the  protrusion  can  be  removed  and  its  effect  in  the  exterior  region  can  be  rep¬ 
resented  by  a  set  of  equivalent  electric  and  magnetic  currents  on  the  surface  of 
the  protrusion.  The  equivalent  electric  and  magnetic  currents  are  related  to  the 
electric  and  magnetic  fields  by 

Js  =  n  X  H,  Ms  =  E  X  n  (1) 

where  n  denotes  the  outward  unit  vector  normal  to  the  surface  of  the  protru¬ 
sion.  Apparently,  the  original  problem  becomes  the  two  equivalent  subproblems, 
depicted  in  Fig.  2.  The  first  is  the  scattering  by  the  large  body  without  the  protru¬ 
sion  and  the  second  is  the  radiation  of  the  equivalent  currents  in  the  presence  of  the 
large  body.  Whereas  the  field  scattered  by  the  large  body  without  the  protrusion 
can  be  calculated  efficiently  and  accurately  using  the  SBR  method  [13]-[15],  the 
calculation  of  the  equivalent  currents  and  their  radiation  is  more  involved.  In  the 
proposed  method,  the  currents  are  calculated  using  the  FEM  in  conjunction  with 
an  absorbing  boundary  condition  (ABC)  and  the  radiation  is  calculated  using  the 
SBR  method  with  the  aid  of  the  reciprocity  theorem.  These  are  discussed  below 
in  more  detail. 
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To  use  the  FEM  with  an  ABC  to  calculate  the  equivalent  electric  and  magnetic 
currents,  we  first  enclose  the  protrusion  in  a  larger  artificial  surface  denoted  as  Fq 
(see  Fig.  3).  We  then  apply  the  second-order  ABC  [16]  to  the  scattered  field  on 


r.: 


d(f>^ 


(2) 


dn  ^  ds^ 

where  (f)=  for  transverse  magnetic  (TM)  or  E^-polarization,  for  trans¬ 

verse  electric  (TE)  or  .ff^j-polarization,  n  denotes  the  outward  unit  vector  normal 
to  Fa,  s  is  the  arc  length  measured  along  Fa,  and 
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where  k{$)  is  the  curvature  of  Fq  at  s.  Since  the  field  scattered  by  the 

protrusion  can  be  considered  as  the  difference  between  the  total  field  {<f>)  and  the 
incident  field  in  the  presence  of  the  large  body  the  ABC  can  be  written  as 


d4>  .  .av 
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where 
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shr 
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This  boundary  condition,  together  with  the  Helmholtz  equation  satisfied  by  (f> 
inside  Fa,  defines  a  unique  boundary- value  problem,  which  can  be  solved  using  the 
FEM  [16].  Once  (f)  inside  Fa  is  found,  the  equivalent  electric  and  magnetic  currents 
on  the  surface  of  the  protrusion  can  be  obtained  using  their  definition  in  (1).  Note 
that  in  (5),  is  calculated  using  the  SBR  method. 

Once  the  equivalent  currents  are  calculated,  their  radiated  field  in  the  presence 
of  the  large  body  can  be  calculated  using  the  SBR  method  with  the  aid  of  the 
reciprocity  theorem  [1],  [2],  [17].  To  be  more  specific,  for  the  TM  polarization,  we 
place  an  infinitely  long  current  filament  Jo  at  the  observation  point  p^.  If  is  far 
from  the  origin,  the  free  space  electric  field  radiated  by  this  current  is  given  by 


E'Jp)  = 

V  SttPo 


(6) 
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where  770  is  the  free-space  intrinsic  impedance.  Apparently,  if  we  choose 


Jo  = - 1 

Vo  ' 


I^JkoPo 

jh 


(7) 


E'(p)  then  becomes  a  plane  wave  incident  from  the  direction  of  the  observation 
point.  If  the  current  filament  is  placed  at  the  observation  point  in  the  presence 
of  the  large  body  without  the  protrusion,  the  field  can  be  calculated  conveniently 
using  the  SBR  method.  Using  the  reciprocity  theorem,  the  field  radiated  by  the 
equivalent  currents  is  then  given  by 


=  1  /  -  Hf^Mt)dl  (8) 

Jo 

where  Fj,  denotes  the  surface  of  the  protrusion,  and  are  the  field  due  to 
the  current  filament  and  are  calculated  using  the  SBR  method.  Similarly,  for  the 
TE  polarization,  the  field  radiated  by  the  equivalent  currents  is  given  by 

Hr(p)  =  (9) 


where  Mo  is  given  by 

Mo  =  (10) 

Note  that  in  the  formulation  described  above,  because  of  the  use  of  an  ABC  on 
Fo,  we  neglected  the  field  scattered  by  the  protrusion,  reflected  and/or  diffracted 
back  to  the  protrusion  by  the  large  body,  and  scattered  by  the  protrusion  again.  In 
most  problems,  this  contribution  is  insigniflcant.  However,  when  the  protrusion  is 
very  close  to  edges  and  reflecting  surfaces,  the  contribution  can  become  significant 
and  its  omission  can  cause  a  substantial  error  in  the  solution.  When  this  happens, 
we  can  use  either  PO  or  the  SBR  method  to  calculate  the  field  radiated  by  the 
equivalent  currents  and  reflected  and/or  diffracted  back  to  the  protrusion.  Using 
this  as  the  secondary  incident  field,  we  can  update  the  equivalent  currents.  This 
process  can  be  repeated  until  there  is  no  significant  change  in  the  values  of  the 
equivalent  currents.  The  field  radiated  by  the  currents  is  then  calculated  using 
either  (8)  or  (9). 
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III.  NUMERICAL  RESULTS 

In  this  section,  we  present  several  examples  to  demonstrate  the  validity  of  the 
hybrid  method  described  above. 

The  first  example  is  a  circular  cylinder  with  a  small  conducting  protrusion. 
The  radius  of  the  circular  cylinder  is  5A  and  the  protrusion  is  lA  wide  and  lA 
high.  The  monostatic  radax  cross  section  (RCS)  is  given  in  Fig.  4.  As  can  be 
seen,  the  hybrid  solution  is  in  good  agreement  with  the  MoM  solution  [18].  In  the 
calculation,  the  absorbing  boundary  is  placed  0.5A  away  from  the  surface  of  the 
protrusion. 

The  second  example  is  a  16Ax8A  rectangular  conducting  cylinder  having  a 
conducting  protrusion  on  its  upper  surface.  The  protrusion  is  0.8A  wide  and  0.8A 
high.  The  monostatic  RCS  is  shown  in  Fig.  5  and  compared  to  the  MoM  solu¬ 
tion.  The  accuracy  here  is  slightly  worse  than  that  in  the  first  example,  and  this 
degradation  is  due  to  the  error  in  the  PO  solution  of  the  scattered  field. 

The  third  example,  shown  in  Fig.  6,  places  a  conducting  protrusion  at  the 
corner  of  the  rectangular  cylinder.  Again,  reasonably  good  agreement  is  obtained. 
Therefore,  the  method  is  not  limited  to  the  protrusions  placed  on  a  locally  flat 
surface. 

The  last  example  is  an  L-shaped  conducting  body  with  a  protrusion  on  the 
surface.  This  problem  differs  from  the  former  ones  in  that  both  incident  and 
scattered  fields  can  have  multiple  bounces.  The  hybrid  solution  obtained  without 
iteration  is  given  in  Fig.  7,  from  which  a  noticeable  error  is  observed.  This  is 
expected  because  the  protrusion  is  close  to  the  reflecting  surface.  The  results 
obtained  using  the  iterative  approach  are  also  shown  in  Fig.  7,  which  demonstrate 
clearly  the  improvement  achieved  by  the  iterative  approach. 
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IV.  CONCLUSION 

In  this  article,  we  presented  the  hybrid  FEM/SBR  method  to  compute  scatter¬ 
ing  by  large  bodies  with  small  protruding  scatterers.  We  first  employed  the  field 
equivalence  principle  to  replace  the  protruding  scatterers  by  a  set  of  equivalent 
electric  and  magnetic  currents.  We  then  used  the  SBR  method  to  calculate  the 
field  scattered  by  the  large  body  without  protrusions  and  the  FEM  to  compute 
the  equivalent  currents.  The  field  radiated  by  the  equivalent  currents  was  also 
calculated  using  the  SBR  method  with  the  aid  of  the  reciprocity  theorem  and  was 
superimposed  to  the  field  scattered  by  the  large  body.  Two-dimensional  problems 
were  given  to  demonstrate  the  feasibility  of  the  proposed  technique.  Its  extension 
to  the  three-dimensional  space  is  straightforward. 
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FIGURE  CAPTIONS 


Fig.  1  Original  problem:  A  large  PEC  body  with  a  small  protrusion. 

Fig.  2  Equivalent  problem:  The  protrusion  is  replaced  by  equivalent  electric  and 
magnetic  currents. 

Fig.  3  Region  for  FEM  calculation:  The  protrusion  is  enclosed  in  an  artificial 
boundary. 

Fig.  4  Comparison  of  the  monostatic  RCS  calculated  by  the  hybrid  SBR/FEM 
and  the  MoM  for  a  circular  cylinder  with  a  conducting  protrusion. 

Fig.  5  Comparison  of  the  monostatic  RCS  calculated  by  the  hybrid  SBR/FEM 
and  the  MoM  for  a  rectangular  cylinder  with  a  conducting  protrusion. 

Fig.  6  Comparison  of  the  monostatic  RCS  calculated  by  the  hybrid  SBR/FEM 
and  the  MoM  for  a  rectangular  cylinder  with  a  conducting  protrusion  at  the 
corner. 

Fig.  7  Comparison  of  the  monostatic  RCS  calculated  by  the  hybrid  SBR/FEM 
and  the  MoM  for  an  L-shaped  cylinder  with  a  conducting  protrusion. 
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Figure  1:  Original  problem:  A  large  PEC  body  with  a  small  protrusion. 


Figure  2:  Equivalent  problem:  The  protrusion  is  replaced  by  equivalent  electric 
and  magnetic  currents. 


Figure  3:  Region  for  FEM  calculation:  The  protrusion  is  enclosed  in  an  artificial 
boundary. 


Theta  (degree) 


Figure  4:  Comparison  of  the  monostatic  RCS  calculated  by  the  hybrid  SBR/FEM 
and  the  MoM  for  a  circular  cylinder  with  a  conducting  protrusion. 


Figure  5:  Comparison  of  the  monostatic  RCS  calculated  by  the  hybrid  SBR/FEM 
and  the  MoM  for  a  rectangular  cylinder  with  a  conducting  protrusion. 


23 


- MOM 

—  Hybrid  Method 
Without  Protrusion 


Figure  7:  Comparison  of  the  monostatic  RCS  calculated  by  the  hybrid  SBR/FEM 
and  the  MoM  for  an  L-shaped  cylinder  with  a  conducting  protrusion. 
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SCATTERING  FROM  A  CYLINDRICALLY  CONFORMAL 
SLOTTED- WAVEGUIDE  ARRAY  ANTENNA 

GuoXin  Fan  and  Jian-Ming  Jin 
Electromagnetics  Laboratory 
Department  of  Electrical  and  Computer  Engineering 
University  of  Illinois  at  Urbana-Champaign 
Urbana,  Illinois  61801-2991 

ABSTRACT 

A  numerical  method  is  developed  to  investigate  electromagnetic  scattering  by  a  cylindrically 
conformal  waveguide-fed  slot  array.  The  problem  is  first  formulated  in  terms  of  integral  equa¬ 
tions  using  the  equivalence  principle.  The  integral  equations  are  then  solved  using  the  method 
of  moments  (MoM)  in  conjunction  with  global  sinusoidal  basis  functions  and  Galerkin’s  testing 
procedure.  The  MoM  solution  requires  the  evaluation  of  the  generalized  admittance  matrices 
involving  various  dyadic  Green’s  functions.  The  slow  convergence  of  the  series  associated  with 
the  summation  of  waveguide  modes  is  accelerated  using  the  Kummer  transformation  and  the 
slow  convergence  of  the  series  associated  with  the  summation  of  the  exterior  modes  is  avoided 
by  using  the  asymptotic  solutions  with  proper  treatment  for  singular  integrals.  The  evaluation 
of  the  excitation  vector  and  scattered  field  is  also  accelerated  using  Watson’s  transformation 
and  asymptotic  solutions.  Numerical  results  are  presented  to  illustrate  the  scattering  character¬ 
istics  of  the  cylindrically  conformal  waveguide-fed  slot  arrays,  such  as  the  effects  of  curvature, 
slot  thickness,  and  waveguide  termination  on  the  radar  cross  section  of  the  arrays. 
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I.  INTRODUCTION 


Slotted-waveguide  array  antennas  are  widely  used  on  modern  airborne  radars  because  of 
their  large  power  handling  capability,  high  efficiency,  light  weight,  compact  structure,  and  low 
side  lobes.  However,  being  highly  efficient  radiating  structures,  by  the  reciprocity  theorem  these 
antennas  are  also  efficient  scatterers,  which  contribute  significantly  to  the  overall  radar  cross 
section  (RCS)  of  the  host  vehicle.  Usual  RCS  reduction  techniques  cannot  be  applied  to  these 
antennas  without  degrading  their  performance.  Bandpass  radomes,  made  of  frequency  selective 
surfaces  (FSS),  can  reflect  nearly  all  the  incident  energy  outside  the  working  frequency  band 
and  allow  the  incident  energy  within  the  working  frequency  band  pass  through  to  reach  the 
antennas.  However,  it  is  at  the  working  frequency  band  that  the  slot  arrays  have  a  significantly 
higher  RCS.  Therefore,  it  is  necessary  to  predict  the  RCS  of  the  slot  arrays  for  applications 
such  as  target  identification,  electromagnetic  compatibility,  and  stealth  technology. 

Because  of  a  large  number  of  slots  and  the  mutual  coupling  between  the  slots  through 
the  waveguides  and  exterior  space,  a  full-wave  analysis  of  slotted  waveguide  arrays  is  very 
difficult.  In  the  past,  most  work  was  focused  on  the  radiation  analysis  of  a  single  slot  [1]-[10], 
one-dimensional  arrays  and  small  two-dimensional  arrays  [11]-[17],  all  in  a  planar  surface.  In 
particular,  Stevenson  [1]  developed  what  is  now  considered  the  classical  theory  for  a  single  slot 
and  Oliner  [2]  presented  a  variational  solution  of  the  problem.  Khac  and  Carson  [3]  employed 
the  method  of  moments  (MoM)  to  seek  a  numerical  solution  to  the  slot  field  using  pulse  basis 
functions  and  point  match  technique.  The  efficiency  of  the  MoM  solution  was  improved  by 
Lyon  and  Sangster  [5]  and  Stern  and  Elliott  [6]  by  using  global  and  piecewise  sinusoidal  basis 
functions  and  Galerkin’s  technique.  The  MoM  was  also  applied  to  tilted  slots  [9],  dielectric- 
covered  slots  [8],  [10],  edge  slots  [18],  and  slots  in  a  sectoral  waveguide  [19].  Recently,  Fan  [20] 
analyzed  cylindrically  conformal  slotted-waveguide  array  antennas  with  a  curved  waveguide  as 
the  feeding  guide  and  sectoral  guides  as  the  radiating  guides.  Whereas  the  literature  for  the 
radiation  analysis  of  the  slot  arrays  is  abundant,  little  work  was  done  on  the  scattering  from 
the  planar  slot  antennas  [21]-[23],  let  alone  the  curved  ones.  Josefsson  [21]  analyzed  scattering 
by  a  single  slot  in  an  infinitely  long  waveguide  and  Chen  and  Jin  [23]  developed  the  MoM 
and  finite  element  method  (FEM)  solutions  for  the  slots  in  the  waveguides  terminated  with 
arbitrary  loads.  No  work  was  found  on  the  analysis  of  the  scattering  characteristics  of  slotted 
waveguide  arrays  on  a  curved  surface.  In  this  paper,  we  present  such  an  analysis. 

The  key  problem  in  the  analysis  of  slotted  waveguide  arrays  is  to  solve  for  the  slot  aperture 
fields.  Among  various  numerical  methods,  the  MoM  is  most  efficient  and  accurate  for  this 
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purpose  because  only  a  small  number  of  basis  functions  are  needed  to  represent  the  aperture 
fields.  However,  the  major  difficulties  in  the  MoM  are  (i)  its  formulation  requires  the  dyadic 
Green’s  functions  for  the  waveguide,  the  exterior  space,  and  the  slot  if  it  has  a  finite  thickness, 
and  (ii)  the  evaluation  of  its  matrix  involves  highly  singular  integrals  and  very  slowly  converging 
series.  In  this  paper,  we  address  all  these  problems.  To  be  more  specific,  we  first  formulate 
the  integral  equations  for  the  problem  and  apply  the  MoM  with  global  basis  functions  and 
Galerkin’s  testing  procedure.  We  then  describe  in  detail  the  evaluation  of  the  MoM  matrix 
involving  various  dyadic  Green’s  functions.  The  slow  convergence  of  the  series  associated  with 
the  summation  of  waveguide  modes  is  accelerated  using  the  Kummer  transformation  and  the 
slow  convergence  of  the  series  associated  with  the  summation  of  the  exterior  modes  is  avoided  by 
using  the  asymptotic  solutions  developed  by  Boersma  and  Lee  [24]  and  Bird  [25].  The  evaluation 
of  the  excitation  vector  and  scattered  field  is  also  accelerated  using  Watson’s  transformation 
and  asymptotic  solutions.  Finally,  we  present  some  numerical  results  to  illustrate  the  scattering 
characteristics  of  the  cylindrically  conformal  slotted  waveguide  arrays,  such  as  the  effects  of 
curvature,  slot  thickness,  waveguide  termination,  and  frequency  on  the  RCS  of  the  arrays. 

I.  THEORY 

In  this  section,  we  describe  in  detail  the  formulation  of  the  problem  and  its  solution  by 
the  MoM.  Particular  attention  is  given  to  the  computation  of  the  elements  of  the  generalized 
admittance  matrices. 

A.  Integral  Equations  and  MoM  Solution 

Consider  the  problem  of  electromagnetic  wave  scattering  by  a  waveguide-fed  slot  array  on  a 
cylindrical  surface  whose  radius  is  ps,  as  depicted  in  Fig.  1.  The  cross  section  of  each  waveguide 
is  an  annular  sector  with  inner  radius  pi,  outer  radius  p2,  and  subtended  angle  A(f>,  and  the 
thickness  of  the  curved  walls  of  the  waveguides  is  t  =  ps-  P2.  Each  waveguide  is  terminated  at 
z  =  zi  and  Z2  with  an  arbitrary  load.  All  radiating  slots  are  longitudinal  slots  cut  in  the  outer 
wall  of  the  waveguides,  having  the  same  width  2w  and  different  length  and  offset  with  respect 
to  the  center-line  of  the  waveguides. 

For  the  ith  slot,  its  inner  and  outer  apertures,  Sf  and  Sf,  divide  the  space  into  three  regions: 
the  waveguide  region  (region  o),  the  region  outside  the  cylinder  (region  6),  and  the  slot  region 
(region  c) ,  as  illustrated  in  Fig.  2.  In  accordance  with  the  equivalence  principle,  the  fields  in  the 
three  regions  can  be  decoupled  by  covering  the  apertures  5/  and  Sf  with  a  perfectly  conducting 
surface  and  introducing  equivalent  magnetic  currents  above  and  below  the  perfectly  conducting 
surface.  Denoting  the  equivalent  magnetic  current  below  Sf  as  Mf  and  that  below  Sf  as  Mf , 
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because  of  the  continuity  of  the  tangential  electric  fields  across  the  apertures,  the  equivalent 
magnetic  current  above  5/  is  -Mf  and  that  above  Sf  is  -Mf .  Therefore,  the  field  in  region  a 
is  due  to  M/,  the  field  in  region  b  is  due  to  — Mf  and  the  incident  field,  and  the  field  in  region 
c  is  due  to  — Mf  and  Mf . 

By  enforcing  the  continuity  of  the  tangential  magnetic  fields  across  Sf  and  Sf,  we  obtain 
the  integral  equations  satisfied  by  the  equivalent  magnetic  currents,  Mf  and  Mf ,  as 


£  h;,.(m{)  +  Hji(Mf)  -  h;,(m?)  =  o 

7 

re  5/ 

(1) 

J 

52  H^.(Mf)  -1-  KUiMf)  -  H=,(Mf)  = 

r€5f 

(2) 
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where  the  subscript  t  denotes  the  tangential  component,  the  summation  in  (1)  is  carried  out 
for  all  inner  apertures  in  the  same  waveguide,  and  the  summation  in  (2)  is  carried  out  for  all 
outer  apertures.  The  magnetic  field  is  related  to  the  surface  magnetic  current  by 

H“(r)  =  (r,  r')  •  M{r')dS'  (3) 

5 

where  a  =  o  for  region  c,  b  for  region  b,  and  cfor  region  c,  respectively,  and  G“(r,  r')  denotes  the 
magnetic-source  magnetic-field  dyadic  Green’s  function  in  the  corresponding  region.  Finally, 
HP’’*  is  the  field  due  to  the  incident  field  in  the  presence  of  the  conducting  cylinder  without 
slots. 

To  seek  a  numerical  solution  of  (1)  and  (2),  we  first  expand  each  equivalent  magnetic  current 
using  the  global  sinusoidal  basis  functions 

Nj’’’ 

=  zMj'^  =  i  £  ^  s'/  (4) 

9=1 


where  a^j  =  qTf/Lj  with  Lj  being  the  length  of  the  yth  slot.  Applying  Galerkin’s  procedure, 
the  integral  equations  can  be  converted  into  the  matrix  equation  given  by 


■[y«((.)  +  y«(c,5';5/)] 

[-y«(c,S/;S/)]  • 

1  1 

(5) 

[y«(6)  +  y«(c,S»;S/)]. 

1 N1 J 

1  UWJ 

where  (o)  j  and  (6)  j  are  the  generalized  admittance  matrix  for  regions  a  and  6,  respec- 

tively.  ^c,  is  the  generalized  admittance  matrix  for  region  c,  in  which  S{  is  the 

aperture  on  which  the  field  point  is  located  and  Sj  is  the  aperture  on  which  the  source  point  is 
located.  Finally,  [Ipi]  is  the  generalized  excitation  vector. 

In  the  remainder  of  this  section,  we  describe  the  evaluation  of  the  elements  of  each  gener¬ 
alized  admittance  matrix  and  excitation  vector,  necessary  for  a  numerical  solution  of  (5). 
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B.  Generalized  Admittance  Matrix  Elements  for  Region  a 

The  admittance  matrix  elements  due  to  the  internal  coupling  in  the  sectoral  waveguide 
(region  a)  are  given  by 


Ypg  (a)  =  JJ  JJ (r,  r')  sin  Opi^  sin  a, dS  dS' 


(6) 
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where  is  the  ii-component  of  the  magnetic-source  magnetic-field  dyadic  Green’s  function 
for  the  sectoral  waveguide,  which  is  derived  in  [26]  and  given  by 
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Also,  in  the  above,  u  =  mir/A4>,  C,p  =  cosv(p,  C^'  =  cosi/^',  /i^n  =  A„i„  are 

the  roots  of  J5[,(ATOnP2)  =  0,  Amn  =  1  -  and  J?<„  and  are  the 

reflection  coeificients  for  mode  (m,  n)  at  2  =  Zi  and  z  =  Z2,  respectively.  For  shorted  terminals 
=  —  1,  and  for  matched  terminals  Rmn  =  0. 

The  integrals  in  (6)  can  be  evaluated  analytically  through  the  transformation  of  variables, 
resulting  in 
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where  P^{m)  =  /^(m)/^(m)  with  /^(O)  =  2w  and 

[“  -h  +  t/)]  m  7^  0 

and  5i  denotes  the  offset  of  the  ith  slot  from  the  center  of  the  waveguide.  Also, 
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^pi  ^mn 

lii{m,n)  =  C{m,n)e^thmn[(zoi-zoA-iLi-Lj)/2] 


1  _  (_i)PeTi/imnl-ij  U  _  (_i)9e=*=^* 


■?oi  <  zoj 


(11) 


(12) 


(13) 


30 


with 


C(m,  n)  = 
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and  ^ot  and  zqj  denote  the  center  of  the  fth  and  jth  slots,  respectively.  Finally,  rJ-^{Tn,n)  can 
be  written  in  three  parts  as 


with 
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While  the  summation  in  (11)  can  be  evaluated  without  any  difficulty  when  i  ^  y,  its 
evaluation  when  i  =  j  involves  a  slowly  converging  series  given  by 


5(m)  =  5]  Bim,  n)I"  =  5i(m)  +  52(m) 
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where 
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To  show  this,  let  us  examine  the  asymptotic  behavior  of  the  series  when  n  — >  00.  Using  the 
asymptotic  expression  of  the  Bessel  functions  for  large  arguments,  we  obtain  the  approximate 
solution  to  the  eigenvalue  equation  Bl{Xmnp2)  =  0  as 
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Clearly,  the  asymptotic  form  of  the  terms  in  5i(m)  is  proportional  to  n~^  and  that  in  52  (m)  is 
proportional  to  n~^.  To  accelerate  the  convergence  of  52 (m),  we  apply  the  Kummer  transform 
[27]  to  (22),  resulting  in 
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The  first  summation  in  (25)  converges  rapidly  and  the  second  summation  can  be  written  in  a 


closed  form  as 
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In  passing,  we  note  that  when  pi  and  />2  oo,  the  expressions  for  sectoral  waveguides 
reduce  to  those  for  the  corresponding  rectangular  waveguides  [26]. 


C.  Generalized  Admittance  Matrix  Elements  for  Region  b 

The  admittance  matrix  elements  due  to  the  external  coupling  in  region  b  are  given  by 

=  r')  sin  api$  sin  agj^'  dS  dS'  (28) 

>  ] 

where  is  the  zz-component  of  the  magnetic-source  magnetic-field  dyadic  Green’s  function 
for  the  conducting  cylinder.  The  rigorous  expressions  of  the  dyadic  Green’s  function  involve 
infinite  series  and  infinite  integrations  of  Hankel’s  functions  [28],  [29].  Although  these  expres¬ 
sions  can  be  converted  into  the  forms  suitable  to  numerical  calculation  [29],  it  is  still  very  time 
consuming  to  compute  the  matrix  elements,  especially  when  the  cylinder  is  large.  Therefore, 
it  is  necessary  to  seek  the  fast  converging  asymptotic  solutions.  Several  different  approximate 
asymptotic  solutions  have  been  developed  in  the  past  [24],  [25],  [30]-[32],  and  each  has  its  ac¬ 
curacy  and  range  of  validity.  A  comparative  study  [29]  shows  that  the  B-L  [24]  and  TSB  [25] 
solutions  offer  the  best  overall  accuracy,  and  these  two  solutions  are  complementary  to  each 
other.  More  specifically,  the  B-L  solution  is  more  accurate  when  the  field  point  is  near  the 
source  or  in  the  paraxial  region  while  the  TSB  solution  is  better  when  the  field  point  is  far 
from  the  source  and  off  axis.  In  this  work,  the  TSB  solution  is  used  for  two  distant  slots  whose 
centers  have  the  same  z-coordinate,  whereas  the  B-L  solution  is  employed  for  all  other  cases. 


32 


When  3^  the  integrals  in  (28)  can  be  evaluated  numerically  without  any  difficulty  [26]. 
However,  when  i  ~  ^zz  a  t  ^  singularity  as  t  ^  0,  and  hence,  a  regularization  is  needed 
to  for  the  evaluation.  To  do  this,  we  first  write  as 


Gl  =  Gt  +  G, 


(29) 


where  (5°^  is  the  iz-component  of  the  dyadic  Green’s  function  for 
plane,  given  by 


an  infinite  conducting  ground 


(30) 


and  Gzz  can  be  considered  as  the  perturbation  due  to  the  curvature  of  the  cylinder  [31],  given 

by 


Gzz  =  Gq (t) I  -  ij  j^cos^ 6  +  q{l-  q){2  —  Z cos^ j 
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where  Go{t)  =  exp(-ifct)/(47rt),  t  =  -  z'}'^  -  {(p  -  6  =  arctan[(z  -  z')/{<f  -  (p% 

q  =  j/kt,  P  =  {kcos*0/2R^)^/H,  and  v  is  the  surface  Fock  function  defined  in  [24]-[26].  As  a 
result,  the  self-admittance  matrix  element  can  be  written  as 

>^i(«’)=i;f(6)+y"(6).  (32) 
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Using  integration  of  part  and  transformation  of  variables,  Fp°“(6)  can  be  written  as  [8] 
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where  t  =  y/u^  -f  and 
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The  integral  in  (33)  becomes  regular  when  evaluated  in  the  polar  coordinates. 

To  evaluate  the  perturbation  term  in  (32),  we  first  let  rj  =  psp,  r)'  =  p^tp',  and  write  it  as 
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Introducing  the  transformations  u  =  ^  v  =  rj  -  r}\  u'  =  ^  +  ^'  -  Li,  and  v'  =  ri  +  rj',  and 
observing  that  Gzz  (w,  u)  is  an  even  function  of  u  and  v,  we  obtain 

=  [l  + ^  i2w- v)Gzziu,v)F{u)dvdu  (36) 
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where 


F(«)  = 
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With  another  transformation  u  =  tcosO,  and  u  =  t  sin  6,  (36)  becomes 
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where  0o  =  tan“^(2ti)/I/,)  and  G«(f,  ^)  has  the  same  form  as  (31).  The  integrals  in  (38) 
become  regular  by  letting  t  =  r^,  which  can  be  evaluated  by  expanding  the  integrand  in  terms 
of  the  power  of  r  and  using  a  Gaussian  quadrature.  The  Fock  function  can  evaluated  using  its 
small-argument  expansions. 


D.  Generalized  Admittance  Matrix  Elements  for  Region  c 

It  can  be  shown  [26]  that  the  magnetic  field  in  the  cavity  (region  c  with  a  sectoral  cross- 
section),  formed  by  the  slot  with  its  two  apertures  covered  by  a  perfectly  conducting  surface, 
due  to  the  surface  magnetic  current  z6{p  —  p')  sin  ap,^  is  given  by 


TT  r  •  (B{p2,p)B{p3,p') 

Hz{p,p)  =  —  sinopi^  < 

jujp  "Hb(p3,/>)B(p3,pO 


p<p' 

p>p' 


where  0^  =  —  al^l  and 


B{p',p)  = 


C=  [B{p2,p')B'{p^,p')  -  B{ps,p')B'{p2,p')]-' 

^  (  Wp>)Yo{l3p)  -  Y>{/3p')JoiM  k  >  api 


mp')KQ{l3p)  -  K'^{^p')h{fip)  k  <  api. 


Using  (39),  we  obtain  the  expressions  for  the  generalized  admittance  matrix  elements  for 
region  c  as 

(«) 

YS(c,S!-.S;)  =  0p,B’(p2,P3) 

^/)  =  teB'(P3,P2) 

(45) 

Note  that  when  p2  and  pz  oo,  the  expressions  above  reduce  to  those  for  a  rectangular  cavity 
[26]. 
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E.  Excitation  Vector  and  Scattered  Field 

Consider  an  incident  plane  wave  of  arbitrary  polarization  and  incidence  angle,  whose  electric 
field  is  given  by 

jjmc  ^  ^  ^  ^inc  gjj,  (45) 

where  ip  is  the  polarization  angle,  denote  the  incidence  angles,  and 

are  the  unit  vectors.  The  primary  field,  in  (2),  is  given  by 


Hr 


P=P3 


=  -Yo  sin  sin  sin  ^ 


(47) 


where 

.r/  \  2  ^  2  r cos n(p 

~  5o  l  +  <^On  iff)  (a:) 

from  which  we  obtain  the  elements  of  the  excitation  vector  as 


(48) 


fl  Hr\  sinapi^dS 
JJ  IP=P3 
Si 


=  4Yo^^  sin  ip  sin  “*^'”7.(0'”")  M(kp3 sin  0*”",  <foi  -  Ay’)  (49) 

pn 


where 


f^^e)=\l-  {-^cos9 


21  ^  J  -  COS  COS  p  =  odd 


sm 


(^cos0) 


p  =  even 


cosnip  .  ,  *  . 

sinc(n  A<p) 


in  which  Ay>  =  tn/ps.  For  narrow  slots,  we  can  use  the  midpoint  integration  to  find 


hi  =  4Yo 


AipLi 
pir 


sin  V’  sin  ‘=°®®‘"Vp(r"=)M(ifcp3sin  (po,-  -  <p’"'). 


(50) 

(51) 

(52) 


When  X  becomes  large,  M{x,ip)  for  the  lit  region  can  be  evaluated  more  efficiently  using  the 
expression  [28] 

M(®,<p)S2e^'*~*'^  (53) 

which  can  be  considered  as  the  geometrical  optics  approximation  and  it  has  the  same  accuracy 
as  Gorianiov’  more  complicated  expression  [33].  For  the  neighborhood  of  the  shadow  boundary 
and  the  shadow  region,  M(x,y>)  can  be  evaluated  using  the  expression  obtained  using  the 
Watson  transformation  and  Fock  theory  [28], 


M(x,(p)~e-M'P-f)(7 


- 1) (I)"' (f-)(l) 


1/31 


(54) 
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where  g{x)  is  the  far-zone  Fock  function,  whose  definition  and  evaluation  are  discussed  in  [26], 
[28],  [34],  and  [35]. 

Once  the  equivalent  magnetic  currents  on  the  slot  apertures  are  found,  we  can  compute  the 
scattered  field  using  the  reciprocity  theorem. 

III.  NUMERICAL  RESULTS 

The  formulation  described  above  has  been  implemented  in  a  computer  program  and  a  variety 
of  numerical  results  have  been  obtained  to  study  the  effects  of  various  factors  on  RCS.  Because 
of  limited  space,  we  present  only  some  typical  results;  more  results  can  be  found  in  [26]. 

For  all  the  examples  considered,  the  cross-section  of  sectoral  waveguides  is  given  by  -f 
P2)A4>  =  21  mm  and  p2  —  Pi  =  10.16  mm.  The  distance  between  the  center-lines  of  the  adjacent 
waveguides  is  22  mm,  measured  on  the  surface  of  the  cylinder.  The  adjacent  slots  on  the  same 
waveguide  are  spaced  Apo/2  apart  in  the  z  direction,  where  A^o  is  the  waveguide  wavelength  at 
the  working  frequency  (the  center  frequency  of  design).  The  shorted  planes  are  placed  Xgo/4 
away  from  the  first  and  last  slots,  respectively.  In  order  to  avoid  the  numerical  overflow  in  the 
computation  of  the  generalized  admittance  elements  for  the  shorted  waveguides  at  the  working 
frequency,  the  position  of  the  shorted  planes  is  displaced  by  Az,  which  can  be  considered  as  an 
acceptable  manufacturing  tolerance.  All  other  parameters,  unless  otherwise  specified,  are  given 
by:  2w  =  1.6  mm,  Li  =  16  mm.  Si  —  1.5  mm.  At  =  0.8  mm,  Az  =  0.01  mm,  V’  =  90°,  and  the 
working  frequency  /o  =  9.1  GHz.  The  normal  at  the  center  of  the  array  is  =  90°  and  =  0° 
in  the  spherical  coordinate  system. 

The  computer  program  was  first  validated  by  comparing  its  solution  to  those  obtained  by 
MoM  and  FEM  [23]  for  a  single  slot  and  for  four  slots  on  a  planar  surface.  Good  agreement 
was  observed  between  the  three  solutions. 

To  demonstrate  the  effect  of  the  curvature  of  the  host  cylinder  on  the  RCS  of  a  conformal 
slot  array.  Figs.  3  and  4  show  the  RCS  for  pi  =  oo,  =  500  cm,  pi  =  200  cm,  pi  =  100  cm,  and 
Pi  =  50  cm  in  both  El-  and  H-planes.  As  can  be  seen,  the  cylinder’s  curvature  has  a  significant 
effect  on  the  patterns  and  values  of  the  RCS  in  the  El-plane.  When  pi  becomes  large,  the  RCS 
approaches  that  for  a  planar  array.  When  pi  decreases,  the  grating  lobes  begin  to  disappear 
and  the  maximum  RCS  decreases.  In  contrast,  the  pattern  of  the  RCS  in  the  H-plane  remains 
the  same  and  only  the  value  of  the  RCS  is  decreased  as  pi  decreases. 

To  further  demonstrate  the  effect  of  the  cylinder’s  curvature  on  the  RCS,  Fig.  5  presents 
the  space  distribution  of  the  RCS  for  pi  =  oo,  pi  =  200  cm,  pi  =  100  cm,  and  pi  =  30  cm.  As 
can  be  seen  clearly,  as  pi  decreases  the  grating  lobes  in  the  El-plane  is  dilated  where  the  pattern 
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in  the  H-plane  remains  the  same. 

To  show  the  eflFect  of  the  waveguide  terminals  on  the  RCS,  Fig.  6  gives  the  RCS  of  a  8  x  8 
conformal  slot  array  {pi  =  100  cm)  when  the  waveguide  terminals  are  matched  and  shorted.  It 
can  be  seen  that  when  the  plane  wave  is  incident  normally  on  the  cylinder  {6  =  90°),  the  RCS 
is  very  small  for  the  array  with  shorted  terminals.  This  is  because  at  the  working  frequency,  the 
wave  entered  the  waveguide  is  reflected  by  the  shorted  terminals  and  the  reflected  wave  cancels 
the  incident  wave  at  the  slot.  Thus,  the  total  electric  field  or  the  total  equivalent  magnetic 
current  is  zero  at  the  aperture  of  the  slot.  An  alternative  explanation  is  that  with  the  shorted 
terminals,  the  equivalent  impedance  looking  at  the  aperture  of  the  slots  is  zero.  Therefore, 
when  the  shorted  planes  are  displaced,  the  equivalent  impedance  will  be  changed  and  so  is  the 
RCS.  This  is  clearly  demonstrated  in  [26]. 

The  elfect  of  the  slot  thickness  on  the  RCS  is  also  investigated  in  [26].  It  is  observed  that, 
at  the  working  frequency,  the  slot  thickness  has  no  effect  on  the  RCS  when  the  waveguides  are 
matched  and  when  the  waveguides  are  shorted,  it  has  a  significant  effect.  This  effect  can  also 
be  explained  using  the  concept  of  the  equivalent  impedance. 

Finally,  Fig.  7  shows  the  frequency  responses  of  a  8  x  8  array  (pi  =  100  cm)  with  matched 
and  shorted  terminals,  respectively.  It  is  very  interesting  to  note  that  the  RCS  of  the  slot  array 
is  almost  identical  no  matter  the  waveguides  are  matched  or  shorted  when  the  frequency  is  not 
close  to  the  working  frequency.  At  the  working  frequency,  the  RCS  is  substantially  different. 
This  implies  that  when  the  frequency  of  the  incident  wave  is  not  close  to  the  working  frequency 
of  the  slot  array,  the  energy  entering  the  waveguides  is  very  trivial  and  the  main  contribution  to 
the  RCS  is  the  wave  scattered  by  the  slots  directly.  Therefore,  in  this  case  the  dominant  factor 
is  the  geometry  of  the  slots,  instead  of  the  structures  behind  the  slots.  Figure  8  shows  the  space 
distribution  of  the  RCS  of  a  16  x  16  array  {pi  =  100  cm)  at  three  different  frequencies. 

IV.  CONCLUSION 

In  this  paper,  a  method  of  moments  (MoM)  solution  was  developed  for  electromagnetic 
scattering  by  a  cylindrically  conformal  waveguide-fed  slot  array.  Numerical  results  were  pre¬ 
sented  to  illustrate  the  scattering  characteristics  of  the  cylindrically  conformal  waveguide-fed 
slot  arrays,  such  as  the  effects  of  curvature,  slot  thickness,  and  waveguide  termination  on  the 
radar  cross  section  of  the  arrays.  It  was  observed  that  (i)  as  the  curvature  of  the  host  cylinder 
increases,  the  grating  lobes  in  the  El-plane  dilate  and  eventually  disappear  and  the  maximum 
RCS  is  reduced;  (ii)  the  internal  structure  of  the  waveguides  such  as  the  waveguide  terminals 
has  negligible  effect  on  the  RCS  when  the  frequency  of  the  incident  wave  is  not  close  to  the 
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working  frequency  of  the  slot  array;  (iii)  at  the  working  frequency  the  RCS  of  the  slot  array  with 
a  shorted  terminals  at  normal  incidence  changes  significantly  with  the  position  of  the  shorted 
terminals;  and  (iv)  at  the  working  frequency  the  slot  thickness  has  no  effect  on  the  RCS  of 
the  array  with  matched  terminals  whereas  for  the  array  with  shorted  terminals  its  effect  is 
noticeable. 
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FIGURE  CAPTIONS 


Figure  1.  Conformal  slotted-waveguide  array. 

Figure  2.  Division  of  three  regions  and  surface  magnetic  currents. 

Figure  3.  The  effect  of  the  cylinder’s  curvature  on  the  RCS  in  the  E)-plane  for  a  16  x  16  slot 

array  with  matched  terminals  at  /  =  9.1  GHz.  (a) - planar, - pi  =  500  cm,  (b)  pi  =  200 

cm,  (c)  Pi  =  100  cm,  (d)  pi  =  50  cm. 

Figure  4.  The  effect  of  the  cylinder’s  curvature  on  the  RCS  in  the  H-plane  for  a  16  x  16  slot 

array  with  matched  terminals  at  /  =  9.1  GHz.  (a) - planar, - pi  =  500  cm,  (b)  pi  =  200 

cm,  (c)  Pi  =  100  cm,  (d)  pi  =  50  cm. 

Figure  5.  RCS  (dBsw)  of  a  16  x  16  slot  array  with  shorted  terminals  at  /  =  9.1  GHz.  (a) 
planar;  (b)  pi  =  200  cm;  (c)  pi  =  100  cm;  (d)  pi  =  30  cm. 

Figure  6.  The  effect  of  the  terminals  on  the  RCS  of  a  8  X  8  conformal  slot  array  (pi  =  100  cm). 
•  •  •  matched  terminals, - shorted  terminals,  (a)  E-plane  pattern;  (b)  H-plane  pattern. 

Figure  7.  RCS  of  a  8  x  8  conformal  slot  array  (pi  =  50  cm)  as  a  function  of  frequency.  0*”^’  =  90°, 

^inc  _  450_ 

Figure  8.  RCS  (dBsw)  of  a  16  x  16  conformal  slot  array  (pi  =  100  cm)  at  three  frequencies, 
(a)  /  =  8  GHz;  (b)  /  =  9.1  GHz;  (c)  /  =  12  GHz. 
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Figure  1:  Conformal  slotted- waveguide  array. 


z 


Figure  2:  Division  of  three  regions  and  surface  magnetic  currents. 


Figure  3:  The  effect  of  the  cylinder’s  curvature  on  the  RCS  in  the  Ei-plane 

for  a  16  X  16  slot  array  with  matched  terminals  at  /  =  9.1  GHz.  (a) - 

planar,  - Pi  =  500  cm,  (b)  pi  =  200  cm,  (c)  pi  =  100  cm,  (d)  pi  =  50 

cm. 


43 


(b) 


(d) 


Figure  4:  The  effect  of  the  cylinder’s  curvature  on  the  RCS  in  the  H-plane 

for  a  16  X  16  slot  array  with  matched  terminals  at  /  =  9.1  GHz.  (a) - 

planar,  -  Pi  =  500  cm,  (b)  pi  =  200  cm,  (c)  pi  =  100  cm,  (d)  pi  =  50 

cm. 
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Figure  5:  RCS  (dBsw)  of  a  16  x  16  slot  array  with  shorted  terminals  at 
/  =  9.1  GHz.  (a)  planar;  (b)  pi  —  200  cm;  (c)  pi  =  100  cm;  (d)  pi  =  30  cm. 
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(b) 

Figure  6:  The  effect  of  the  terminals  on  the  RCS  of  a  8  x  8  conformal  slot 

array  {pi  =  100  cm).  •  •  •  matched  terminals,  -  shorted  terminals,  (a) 

E-plane  pattern;  (b)  H-plane  pattern. 
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Figure  7:  RCS  of  a  8  x  8  conformal  slot  array  {pi  =  50  cm)  as  a  function  of 
frequency,  (a)  =  90°,  =  0°;  (b)  =  90°,  =  45°. 
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APPENDIX  3 


Hybrid  MoM/SBR  Method  to  Compute 
Scattering  from  a  Slot  Array  Antenna 
in  a  Complex  Geometry 

Andrew  D.  Greenwood  and  Jian-Ming  Jin 
Center  for  Computational  Electromagnetics 
University  of  Illinois 
Urbana,  Illinois  61801 


1  Introduction 

The  presence  of  a  slotted  waveguide  array  antenna  on  a  radar  target  may  have 
a  significant  contribution  to  the  overall  radar  cross-section  (RCS)  of  the  target. 
Therefore,  the  computation  of  the  RCS  should  include  the  scattering  from  the 
slot  array.  Recently,  a  method  of  moments  (MoM)  procedure  has  been  introduced 
to  compute  the  scattering  from  a  cylindrically  conformal  slotted-waveguide  array 
antenna  [1,2].  However,  this  procedure  does  not  take  into  account  the  geometry 
in  which  the  slot  array  is  located.  If  the  slot  array  is  located  in  a  complex,  three 
dimensional  (3-D)  geometry,  the  MoM  cannot  efficiently  account  for  the  effect  of 
the  geometry.  A  more  efficient  method  to  compute  the  scattering  from  a  large,  3-D 
body  is  the  high  frequency  shooting  and  bouncing  ray  (SBR)  method.  However, 
this  method  cannot  accurately  account  for  the  slots,  each  of  which  is  typically 
smaller  than  an  electromagnetic  wavelength  in  size.  In  this  paper,  the  MoM  com¬ 
putation  of  the  scattering  from  a  slot  array  is  hybridized  with  the  SBR  method  to 
compute  the  electromagnetic  scattering  from  a  large,  3-D  target  which  includes  a 
slot  array  antenna. 

The  basis  of  the  hybrid  method  is  the  field  equivalence  principle,  which  allows 
the  scattering  geometry  to  be  decomposed  into  separate  regions.  The  MoM  is 
applied  to  the  slotted  waveguides,  while  the  SBR  method  is  applied  to  the  region 
outside  the  waveguides,  which  includes  the  complex,  3-D  target.  By  using  the 
hybrid  method,  the  scattering  from  a  large,  3-D  target,  which  includes  a  slotted- 
waveguide  array  antenna,  can  be  efficiently  and  accurately  computed. 
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(a)  Large  target 


Region  I 
-M  ^  (outside) 


Region  HI 
(waveguide) 


Region  n 
(slot) 


(b)  Local  slot  geometry 

Figure  1:  Example  of  a  complex,  3-D  target  with  a  slot  array  antenna  and  the 
local  slot  geometry. 

The  remainder  of  this  paper  is  divided  into  five  sections.  Section  2  describes 
the  formulation  of  the  problem,  including  the  use  of  the  MoM,  the  use  of  the  SBR 
method,  and  techniques  to  decouple  the  computations  of  the  two  methods.  Section 

3  describes  briefly  how  the  method  has  been  tested,  and  Section  4  gives  some 
numerical  results  which  show  the  capability  of  the  method.  The  results  in  Section 

4  also  demonstrate  the  need  to  include  the  slot  array  in  scattering  computations. 
Finally,  Section  5  gives  a  brief  conclusion. 

2  Formulation 

Consider  the  example  target  shown  in  Figure  la.  The  target  is  complex  and  3- 
D,  and  it  includes  a  slotted  waveguide  array  antenna  on  its  surface.  The  slotted 
waveguide  array  antenna  may  be  planar,  or  it  may  conform  to  the  surface  of  a 
cylinder.  The  first  step  to  compute  the  scattering  from  this  target  is  to  analyze 
the  slotted  waveguides  using  the  MoM.  Then,  the  scattering  from  the  target  with 
the  slot  apertures  covered  by  perfect  electric  conductor  (PEC)  is  computed  using 
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the  SBR  method.  During  the  SBR  calculation,  the  incident  field  on  the  slot  array 
antenna  is  computed  and  stored.  This  incident  field  is  combined  with  the  MoM 
analysis  to  find  an  equivalent  magnetic  current  on  the  outer  aperture  of  each  slot. 
Finally,  the  radiation  of  these  equivalent  magnetic  currents  in  the  presence  of  the 
complex,  3-D  target  is  computed  using  the  reciprocity  theorem.  This  result  is 
added  to  the  previously  computed  SBR  scattering  result. 


2.1  Use  of  MoM 

The  first  step  in  the  formulation  of  the  problem  is  to  analyze  the  slotted  waveguides 
using  the  MoM.  There  are  two  main  steps  in  the  application  of  the  MoM.  First, 
the  problem  must  be  described  in  terms  of  an  integral  equation.  Then,  the  integral 
equation  is  discretized  to  find  a  numerical  solution.  The  steps  are  outlined  here, 
and  more  detail  is  given  in  [1,2]. 

To  derive  the  integral  equation,  the  apertures  of  each  slot  are  first  covered 
with  PEC,  and  equivalent  magnetic  currents  over  each  aperture  are  introduced. 
Figure  lb  depicts  the  situation  for  the  slot.  The  region  outside  of  the  antenna 
is  denoted  Region  I,  the  region  inside  the  slot  is  Region  II,  and  the  region  outside 
of  the  slot  but  inside  the  waveguide  is  Region  III.  An  equivalent  magnetic  current 
M“  is  introduced  on  the  inside  of  the  outer  slot  aperture  (between  Regions  I  and 
II),  and  the  equivalent  current  is  introduced  on  the  waveguide  side  of  the 
inner  aperture  (between  Regions  II  and  III).  Because  the  electric  field  must  be 
continuous  across  each  aperture,  -M?  must  be  introduced  on  the  outside  of  the 
outer  aperture,  and  — must  be  introduced  on  the  slot  side  of  the  inner  aperture. 
Note  that  when  the  analysis  is  complete,  — M“  are  the  currents  that  radiate  in  the 
presence  of  the  complex,  3-D  body  as  discussed  above. 

To  derive  the  integral  equation,  the  continuity  of  the  tangential  magnetic  field 
across  each  aperture  is  enforced.  Denoting  the  tangential  magnetic  field  in  Region 
III  on  the  slot  aperture  due  to  the  magnetic  current  on  the  aperture  as 
H™(Mj),  the  following  must  hold  on  each  inner  aperture: 

Y:  Hj^'(M5)  -K  -  H^^(M“)  =  0.  (1) 

3 

Further,  denoting  the  tangential  incident  field  on  the  slot  aperture  as  H®f^, 

E  +  HS(Mf )  -  H”  (M*)  =  H??"  (2) 
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must  hold  on  each  outer  slot  aperture.  Note  that  the  incident  fields  are  calculated 
using  the  SBR  method,  and  the  magnetic  field  due  to  a  magnetic  current  is  found 
from 

H"(M)  =  jj  W{r,  r')  •  M{r')dS'  (3) 

5 

where  a  is  I,  II,  or  III,  depending  on  the  region  of  interest,  G  (r,  r')  is  the  magnetic 
source-magnetic  field  dyadic  Green’s  function  in  the  appropriate  region,  and  r 
corresponds  to  the  point  at  which  the  magnetic  field  is  to  be  evaluated.  Combining 
Equations  1,  2,  and  3  gives  an  integral  equation  for  the  magnetic  currents. 

The  second  main  step  in  application  of  the  MoM  is  to  discretize  the  integral 
equation  to  find  a  numerical  solution  for  the  currents.  To  accomplish  this  step,  the 
currents  are  expanded  in  terms  of  sinusoidal  basis  functions.  Defining  ^  to  be  the 
direction  parallel  to  the  lengths  of  the  slots  and  using  a  local  coordinate  system 
in  which  =  0  at  one  end  of  the  slot,  the  current  on  the  slot  aperture  is 
expanded  as 

(4) 

where  N  is  the  number  of  terms  in  the  expansion,  and  (3  represents  a  for  the 
current  on  the  outer  aperture  or  h  for  the  current  on  the  inner  aperture.  Equation 
4  is  valid  for  points  on  the  slot  aperture;  for  points  outside  of  the  aperture,  the 
expansion  is  defined  to  be  zero.  Assuming  the  width  of  a  slot  is  much  less  than  its 
length,  the  ^  component  of  the  current  is  the  only  component  of  interest. 

Substituting  the  expansion  given  in  Equation  4  into  the  integral  equation  allows 
the  integral  equation  to  be  converted  to  a  matrix  equation  which  can  be  solved 
numerically.  For  the  more  details  about  solving  the  integral  equation,  the  reader 
is  referred  to  [1,2].  However,  one  important  step  that  should  be  mentioned  here  is 
the  derivation  of  the  dyadic  Green’s  functions  for  the  various  regions.  The  Green’s 
functions  given  in  [1, 2]  for  Regions  II  and  III  are  applicable  to  the  present  problem. 
For  Region  I,  the  dyadic  Green’s  function  can  be  written  as 

V\v,  t')  =  r')  -h  r').  (5) 

The  Green’s  function  given  in  [1, 2]  for  the  exterior  region  corresponds  to  G  ^  (r,  r'), 
=difr 

and  G  (r,  r')  is  a  perturbation  term  due  to  diffraction  and  reflection  by  the 
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■““diff 

complex  target  in  which  the  slot  array  is  embedded.  Neglecting  G  (r,  r')  neglects 
fields  which  are  scattered  by  the  slots,  diffracted  or  reflected  by  the  large  body 
back  to  the  slots,  and  scattered  by  the  slots  again  [3].  These  fields  are  usually  an 
insignificant  part  of  the  scattering,  and  this  term  is  neglected  in  the  computations. 
Thus,  the  Green’s  function  given  in  [1,2]  for  Region  I  is  used  for  the  present 
problem. 

2.2  Use  of  SBR 

As  previously  mentioned,  the  MoM  is  used  to  analyze  the  slot  array  antenna  while 
the  SBR  method  is  used  for  the  remainder  of  the  problem.  Thus,  there  are  three 
main  tasks  to  be  accomplished  by  the  SBR  method;  to  compute  the  scattering 
from  the  complex,  3-D  target,  to  compute  the  incident  magnetic  fields  on  the  slot 
apertures,  and  to  compute  the  radiation  of  the  equivalent  currents  on  the  slot 
apertures  in  the  presence  of  the  complex,  3-D  target.  In  all  of  these  cases,  the  slot 
apertures  are  covered  with  PEC. 

In  the  SBR  method,  a  dense  grid  of  rays,  corresponding  to  a  plane  wave,  is 
launched  toward  the  target,  and  each  ray  is  traced  as  it  bounces  around  the  target. 
The  bounces  are  governed  by  Geometrical  Optics  (GO),  and  as  each  ray  leaves  the 
target,  its  contribution  to  the  scattering  is  computed  by  a  Physical  Optics  (PO) 
integration.  If  more  accuracy  is  desired,  the  first  order  edge  diffracted  terms  are 
computed  using  the  Geometrical  Theory  of  Diffraction  (GTD)  and  added  to  the 
result  [3-6].  For  the  present  problem,  this  SBR  procedure  is  followed  to  compute 
the  scattering  from  the  complex,  3-D  target  with  the  slot  apertures  closed  by  PEC. 
The  SBR  procedure  is  implemented  using  the  XPATCH  software  package  [4,5]. 

The  incident  magnetic  field  on  the  slot  apertures  is  computed  using  SBR  at  the 
same  time  the  scattering  from  the  complex,  3-D  target  is  computed.  While  tracing 
the  rays  to  find  the  scattering,  some  rays  will  hit  on  or  near  a  slot  aperture.  The 
field  contributions  from  each  of  these  rays  are  combined  with  appropriate  phase 
shifts  to  find  the  incident  magnetic  field  on  each  slot  aperture.  The  incident  mag¬ 
netic  fields  on  the  slot  apertures  are  used  by  the  MoM  to  compute  the  equivalent 
magnetic  currents  on  the  apertures. 

The  remaining  step  in  the  problem  is  to  compute  the  radiation  of  the  magnetic 
currents  in  the  presence  of  the  large  body.  The  SBR  method  together  with  the 
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reciprocity  theorem  is  employed  for  this  task  [3, 6].  Consider  an  infinitesimal  dipole 
placed  at  the  scattering  observation  point.  If  the  target  containing  the  slot  array 
is  in  the  far  field  of  the  dipole,  the  dipole  launches  a  plane  wave  toward  this 
target.  Recall  that  for  the  SBR  method,  the  grid  of  rays  launched  toward  the 
target  corresponds  to  a  plane  wave.  Note  also  that  the  reciprocity  theorem  states 

j IJ  E^'°‘  ‘JdV  =  Jj  ^  •  M“  (6) 

V  s 

where  is  the  incident  field  on  the  slot  apertures  due  to  the  dipole  at  the 

scattering  observation  point,  M“  is  the  current  on  the  outer  slot  apertures,  which 
is  found  using  the  MoM,  is  the  radiation  due  to  — M“,  and  J  is  the  dipole 
current.  Thus,  if  the  dipole  current  (J)  is  appropriately  chosen  and  monostatic 
scattering  is  being  computed,  all  components  to  find  using  reciprocity  are 
computed  already.  If  bistatic  scattering  results  are  desired,  resulting  from  a 

dipole  at  the  scattering  observation  point  must  be  computed  first,  then  can 
be  computed. 

2.3  Decoupling  the  MoM  from  the  SBR  Method 

As  they  are  presented  in  Section  2.1,  the  MoM  computations  are  coupled  to  the 
SBR  method  computations.  This  is  due  to  the  fact  that  the  incident  magnetic  field 
on  the  slot  apertures,  which  is  computed  using  the  SBR  method,  is  required  for  the 
MoM  computations.  To  avoid  having  to  repeat  the  MoM  computations  in  order 
to  analyze  the  scattering  from  many  different  incidence  angles,  it  is  desirable  to 
decouple  the  computations  of  the  two  methods.  There  are  two  ways  of  doing  this. 
The  first  method  preserves  the  coupling  interactions  between  different  slots;  the 
second  involves  an  approximation  which  neglects  the  coupling  between  different 
slots  to  achieve  lower  computational  complexity. 

To  decouple  the  MoM  computations  from  the  SBR  computations  while  preserv¬ 
ing  the  coupling  between  the  various  slots,  the  incident  magnetic  field  on  the  slot 
apertures  can  be  expanded  in  terms  of  basis  functions.  Assuming  that  the  width 
of  a  slot  is  much  less  than  its  length,  the  component  of  the  incident  magnetic  field 
along  the  length  of  a  slot  is  the  only  component  of  interest.  A  convenient  basis 
set  is  the  set  of  pulse  ba^is  functions,  where  the  basis  function  is  defined  to  be 
one  on  the  slot  aperture  and  zero  elsewhere.  There  will  be  n  basis  functions 
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in  the  set,  where  n  is  the  number  of  slots  in  the  array.  The  magnetic  currents  on 
each  slot  aperture  are  then  computed  with  the  incident  field  on  the  slot  array  set 
equal  to  each  of  the  n  basis  functions  in  turn.  A  matrix-vector  multiply  is  then 
carried  out  during  the  SBR  computations.  This  matrix-vector  multiply  converts 
the  incident  magnetic  fields  on  the  slot  apertures  to  the  equivalent  currents  on  the 
apertures. 

The  second  method  of  decoupling  the  MoM  computations  from  the  SBR  com¬ 
putations  neglects  the  coupling  between  the  individual  slots.  One  slot  on  the  slot 
array  is  chosen,  and  it  is  assumed  that  this  slot  is  the  only  slot  present.  The  MoM 
computation  is  carried  out  with  a  magnetic  field  of  unit  amplitude  on  this  slot,  and 
the  result  is  a  magnetic  current  for  the  slot.  It  is  then  assumed  that  all  of  the  slots 
in  the  array  are  equivalent;  the  magnetic  current  on  each  slot  is  set  equal  to  the 
magnetic  field  on  that  slot  times  the  single  magnetic  current  which  is  computed  by 
the  MoM.  This  approximation  reduces  the  computational  complexity  and  the  sizes 
of  data  files.  However,  it  does  not  produce  accurate  results  when  the  frequency  is 
near  the  working  frequency  of  the  slot  array.  This  is  demonstrated  in  Section  4. 

3  Testing 

Before  using  any  new  numerical  technique,  the  technique  should  be  tested  against 
existing  techniques  to  ensure  its  validity.  The  validity  of  the  MoM  computation 
involving  the  coupling  between  the  different  slots  in  the  array  is  validated  by 
comparison  with  previous  MoM  and  finite  element  method  (FEM)  techniques  [2, 7]. 
The  SBR  method  is  also  validated  through  extensive,  previous  testing  [4,5].  The 
hybrid  technique  is  validated  by  comparison  with  a  previous  hybrid  method  to 
compute  the  scattering  from  complex  targets  with  cracks  and  cavities  on  their 
surfaces  [3].  This  is  accomplished  by  considering  a  waveguide  with  a  single  slot 
on  its  surface.  The  slotted  waveguide  is  placed  an  a  large  plate,  and  the  problem 
is  modeled  both  with  the  hybrid  MoM/SBR  method  discussed  in  this  paper  and 
with  the  hybrid  FEM/SBR  method  presented  in  [3].  The  two  solutions  show  good 
agreement. 
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Figure  2:  Configuration  of  the  slots  on  the  surface  of  a  waveguide. 

4  Numerical  Examples 

To  show  the  validity  and  utility  of  the  proposed  technique,  several  numerical  re¬ 
sults  are  presented.  For  all  of  the  numerical  examples,  the  slot  array  contains  16 
waveguides  with  16  slots  on  each  waveguide,  and  the  array  is  designed  to  radiate 
at  9.1  GHz.  In  addition,  the  following  parameters  apply  to  all  of  the  examples 
presented;  the  upper  waveguide  wall  in  which  the  slots  are  cut  is  0.08  cm  thick, 
the  waveguides  are  separated  by  walls  0.1cm  thick,  each  slot  is  1.6  cm  long  and 
0.16  cm  wide,  and  the  slots  are  positioned  on  the  waveguide  surface  as  shown  in 
Figure  2,  where  the  offset  of  each  slot  from  the  center  of  the  waveguide  is  0.15 
cm.  Unless  otherwise  noted,  the  coupling  between  individual  slots  in  the  array  is 
included  in  the  results. 

The  first  example  is  a  planar  slot  array  which  is  in  a  simple  ground  plane 
geometry.  The  waveguides  are  2.230  cm  wide  by  1.016  cm  high,  the  slot  centers 
are  2.444  cm  apart,  and  the  first  and  last  slots  centers  are  1.222  cm  from  the  ends 
of  the  waveguides.  Thus,  the  entire  slot  array  and  the  ground  plate  are  37.3  cm 
wide  by  39.1  cm  long.  In  Figure  3,  the  radar  cross-section  (RCS)  of  the  plate  with 
the  slots  is  superimposed  on  the  scattering  from  the  plate  alone.  The  scattering 
frequency  is  9.1  GHz,  which  is  the  working  frequency  of  the  slot  array.  Figure  3 
shows  results  in  both  the  -plane  and  the  £-plane  and  for  waveguides  which  are 
terminated  both  with  matched  loads  and  with  short  circuits.  For  some  incidence 
angles,  the  slot  array  has  a  dominant  effect  on  the  scattering. 

The  second  example  is  a  slot  array  on  a  cylinder  with  a  nose  cone.  The  radius 
of  the  cylinder  is  16.096  cm,  and  the  length  without  the  nose  cone  is  100  cm.  The 
nose  cone  is  30  cm  long.  The  waveguide  cross-sections  are  sectoral  in  shape  and 
are  1.016  cm  thick.  Along  the  slotted  surface,  the  waveguides  are  2.230  cm  wide. 
The  slots  are  2.573  cm  apart,  and  the  first  and  last  slots  are  1.287  cm  from  the 
ends  of  the  waveguides.  The  entire  slot  array  is  37.3  cm  along  the  circumference  of 
the  cylinder  and  41.2  cm  along  the  axis  of  the  cylinder.  In  Figure  4,  the  i?-plane 
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(a)  if -plane,  matched  waveguide 
loads 


(c)  /T-plane,  short-circuit  waveguide 
loads 


(b)  JS-plane,  matched  waveguide 
loads 


(d)  £'-plane,  short-circuit  waveguide 
loads 


Figure  3:  RCS  of  a  planar  slot  array  on  a  ground  plate  at  9.1  GHz,  the  working 
frequency  of  the  slot  array. 

RCS  of  the  cylinder  alone  and  the  RCS  of  the  cylinder  with  the  slot  array  are 
compared.  Again,  the  scattering  frequency  is  9.1  GHz,  the  working  frequency  of 
the  slot  array,  and  again,  there  are  scattering  directions  for  which  the  slot  array 
dominates  the  return. 

The  next  example  is  intended  to  show  the  effect  of  the  uncoupled  slot  approx¬ 
imation  which  was  discussed  in  Section  2.3.  Figure  5  shows  the  RCS  of  the  same 
geometry  considered  in  the  second  example,  but  as  a  function  of  frequency.  The 
incident  direction  is  40°  in  the  /f-plane.  The  RCS  computed  considering  the  cou¬ 
pling  between  individual  slots  is  plotted  with  the  RCS  computed  by  neglecting  the 
slot  coupling.  The  approximation  neglecting  slot  coupling  is  reasonably  accurate 
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Elevation  Angle  (degrees)  Elevation  Angle  (degrees) 


(a)  Matched  waveguide  loads  (b)  Short-circuit  waveguide  loads 

Figure  4:  RCS  of  a  conformal  slot  array  on  a  cylinder  with  a  nose  cone  at  9.1  GHz, 
the  working  frequency  of  the  slot  array. 


Frequency  (GHz)  Frequency  (GHz) 


(a)  Matched  waveguide  loads  (b)  Short-circuit  waveguide  loads 

Figure  5:  RCS  of  a  conformal  slot  array  on  a  cylinder  with  a  nose  cone.  The 
scattering  is  computed  with  and  without  including  the  coupling  between  individual 
slots.  Near  the  working  frequency  of  the  slot  array  (9.1  GHz),  the  slot  coupling 
must  be  included. 
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away  from  the  working  frequency  of  the  slot  array  antenna,  but  there  is  significant 
error  near  the  working  frequency.  Thus,  this  approximation  must  be  applied  with 
care. 

The  final  example  shows  the  usefulness  of  the  method.  The  slot  array  antenna 
from  the  first  example  is  mounted  on  the  belly  of  an  f309  aircraft,  with  the  lengths 
of  the  slots  perpendicular  the  the  length  of  the  aircraft  body  (see  Figure  6a).  Figure 
6b  shows  the  W-polarized  range  profile  of  the  airplane  both  with  and  without  the 
slot  array.  The  range  profile  is  the  time  domain  response  to  an  incident  sine  pulse. 
The  sine  pulse  in  this  example  has  a  center  frequency  of  10  GHz  and  a  bandwidth 
of  2  GHz,  and  the  slot  array  has  matched  waveguide  loads.  The  slot  scattering 
dominates  the  range  profile. 

5  Conclusion 

A  hybrid  MoM/SBR  method  is  developed  to  compute  the  scattering  from  a  com¬ 
plex,  3-D  target  with  a  slotted  waveguide  array  antenna.  Because  the  target  is 
large  and  3-D,  the  MoM  alone  cannot  efficiently  compute  the  scattering,  and  be¬ 
cause  the  slots  on  the  waveguides  are  small  features,  the  SBR  method  alone  is 
not  accurate.  The  hybrid  method  combines  the  two  individual  methods  in  such  a 
manner  that  the  scattering  can  be  efficiently  and  accurately  computed.  In  the  hy¬ 
brid  method,  the  MoM  is  used  to  model  the  details  of  the  slot  array,  and  the  SBR 
method  is  used  to  model  the  electromagnetic  interactions  with  the  large,  complex 
target.  The  method  is  validated  by  comparison  to  previously  published  methods. 
Numerical  examples  show  the  need  to  include  a  slot  array  model  when  computing 
the  scattering  from  a  complex  target  with  a  slotted  waveguide  array.  The  examples 
also  illustrate  the  capability  of  the  method. 
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(c)  Range  profile  with  slot  array. 

Figure  6:  Range  profile  of  an  f309,  W-polarization,  10  GHz  center  frequency,  2 
GHz  bandwidth.  The  slot  array  has  matched  waveguide  loads. 
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APPENDIX  4 


On  the  Formulation  of  Hybrid  Finite-Element  and 
Boundary-Integral  Method  for  3D  Scattering 
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Abstract-This  paper  studies  in  detail  a  variety  of  formulations  for  the  hybrid 
finite-element  and  boundary-integral  (FE-BI)  method  for  three-dimensional  (3D) 
electromagnetic  scattering  by  inhomogeneous  objects.  It  is  shown  that  the  effi¬ 
ciency  and  accuracy  of  the  FE-BI  method  depend  highly  on  the  formulation  and 
discretization  of  the  boundary-integral  equation  (BIE)  used.  A  simple  analysis  of 
matrix  condition  identifies  the  efficiency  of  the  different  FE-BI  formulations  and 
an  analysis  of  weighting  functions  shows  that  the  traditional  FE-BI  formulations 
cannot  produce  accurate  solutions.  A  new  formulation  is  then  proposed  and  nu¬ 
merical  results  show  that  the  resulting  solution  has  a  good  efficiency  and  accuracy 
and  is  completely  immune  to  the  problem  of  interior  resonance. 
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I.  INTRODUCTION 

The  hybrid  finite-element  and  boundary-integral  (FE-BI)  method  is  a  powerful 
numerical  technique  for  computing  scattering  by  inhomogeneous  objects.  The 
method  first  divides  the  problem  into  an  interior  and  exterior  problems.  The  field 
in  the  interior  region  is  formulated  using  the  finite-element  method  (FEM),  and  the 
field  in  the  exterior  region  is  represented  by  a  boundary-integral  equation  (BIE). 
The  interior  and  exterior  fields  are  then  coupled  by  the  field  continuity  conditions. 

The  hybrid  FE-BI  method  has  been  first  applied  to  two-dimensional  (2D)  scat¬ 
tering  problems  [l]-[6]  and  later  extended  to  more  challenging  three-dimensional 
(3D)  scattering  problems  [7]-[14].  To  be  more  specific,  Paulsen  et  al  [7]  developed 
the  first  FEi-BI  formulation  for  a  general  3D  scattering  problem,  which  employed 
node-based  FEM  to  discretize  the  interior  fields  and  used  either  the  electric-field 
integral  equation  (EFIE)  or  the  magnetic-field  integral  equation  (MFIE)  as  BIE 
to  represent  the  exterior  field.  The  formulation,  however,  exhibited  two  major 
drawbacks.  First,  it  inherited  all  the  difficulties  caused  by  the  use  of  node-based 
elements  to  discretize  the  electric  and  magnetic  fields  directly  [15].  These  diffi¬ 
culties  include  the  treatment  of  dielectric  interfaces  and  sharp  conducting  edges 
and  corners  and  the  appearance  of  spurious  solutions.  Second,  it  failed  at  the 
interior  resonant  frequencies,  which  are  defined  as  the  resonant  frequencies  of  a 
cavity  formed  by  covering  the  surface  where  BIE  applies  with  a  perfect  conductor 
and  filling  its  interior  with  the  exterior  medium.  The  first  difficulty  was  removed 
by  the  use  of  edge-based  FEM  [8],  [9],  [12]-[14]  and  the  second  difficulty  was  alle¬ 
viated  by  the  use  of  the  combined  field  integral  equation  (CFIE),  which  is  a  linear 
combination  of  EFIE  and  MFIE  [10]-[12],  [14]. 

Although  the  FE-BI  method  with  the  implementation  of  edge-based  elements 
and  CFIE  is  remarkably  more  powerful  than  other  numerical  techniques  in  dealing 
with  inhomogeneous  objects,  it  still  has  a  bottleneck  which  is  the  full  matrix  gen¬ 
erated  by  BIE.  As  pointed  out  in  [16],  this  bottleneck  severely  limits  the  capability 
of  the  FEi-BI  method  in  dealing  with  large  objects.  Although  this  problem  can  be 
circumvented  in  some  special  problems  [9],  [16]  or  paxtially  alleviated  using  special 
surfaces  to  separate  the  interior  and  exterior  regions  [11],  [14],  no  efficient  method 
has  been  developed  for  general  3D  problems  so  far. 

Our  renewed  interest  in  the  FE-BI  method  originated  from  the  recent  devel¬ 
opment  of  the  fast  multipole  method  (FMM)  [17]  and  the  multilevel  fast  multipole 
algorithm  (MLFMA)  [18].  Our  objective  is  to  apply  MLFMA  to  BIE  to  completely 
remove  the  bottleneck  in  the  FE-BI  method  for  general  3D  problems.  During  the 
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course  of  pursuing  this  goal,  we  have  encountered  several  problems  associated  with 
the  efficiency  and  accuracy  of  the  FE-BI  method  implemented  using  the  edge-based 
elements  and  CFIE.  This  paper  reports  our  study  of  these  problems. 

In  this  paper,  we  first  formulate  the  general  FE-BI  method  for  3D  scattering 
problems.  We  then  show  that  there  are  several  different  approaches  to  the  dis¬ 
cretization  of  CFIE,  yielding  solutions  of  different  efficiency  and  accuracy.  How¬ 
ever,  none  of  the  traditional  approaches  produces  satisfactory  results.  The  cause  is 
determined  and  a  new  formulation  is  proposed.  Furthermore,  we  show  that,  con¬ 
trary  to  the  common  belief,  not  all  CFIE  formulations  are  immune  to  the  problem 
of  interior  resonance;  however,  the  new  one  is. 


II.  FORMULATION  AND  ANALYSIS 

Consider  the  problem  of  electromagnetic  wave  scattering  by  an  arbitrarily- 
shaped  and  inhomogeneous  body  characterized  by  relative  permittivity  and  per¬ 
meability  (e,.,  Hr),  which  can  be  complex  if  the  body  is  lossy.  To  solve  this  problem 
using  the  FE-BI  method,  we  first  introduce  an  artificial  surface  S  (which  can  be 
the  surface  of  the  body)  to  enclose  the  body  and  divide  the  problem  into  an  in¬ 
terior  and  exterior  ones.  The  field  inside  S  can  be  formulated  into  an  equivalent 
variational  problem  with  the  functional  given  by  [15] 


^(E)  =  r 

^  ^  2jv 


■(V  X  E)  •  (V  X  E)  -  fc^e^E  •  E 


-t-jfcoJ^(E  xH).nd5 


dV 


(1) 


where  V  denotes  the  volume  enclosed  by  S,  n  denotes  the  outward  unit  normal 
to  S,  ko  is  the  free-space  wavenumber,  and  H  =  ZqH  with  Zq  being  the  free- 
space  intrinsic  impedance.  Using  FEM  with  edge  elements,  we  obtain  the  matrix 
equation 


Kii  Kis  0 
Ksi  Kss  B 


(2) 


where  {Ej}  is  a  vector  containing  the  discrete  electric  fields  inside  V,  {E5}  and 
{JT5}  are  the  vectors  containing  the  discrete  electric  and  magnetic  fields  on  S, 
respectively.  Furthermore,  [Kn],  [A'/s],  [Ksi],  [lifss],  and  [B]  are  sparse  matrices 
and,  in  particular,  [Kji]  and  [/i'ss]  are  symmetric  and  [/if/s]  =  [KsiV,  where  the 
superscript  T  denotes  the  transpose. 
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Equation  (2)  cannot  be  solved  unless  a  relation  between  {£^5}  and  {^’5}  is 
specified.  Such  a  relation  is  provided  by  BIE  for  the  exterior  field,  whose  dis¬ 
cretization  yields 

(PKEs}  +  (0]{5s}  =  {i.}  (3) 

where  {6}  is  a  vector  related  to  the  incident  field.  Combining  (2)  and  (3),  we 
obtain  the  complete  system 

'  Kn  Kis  OU  El  \  f  0  ’ 

Ksi  Kss  B  <  Es  ^  ^  0  ^  (4) 

.0  P  q\[  Hs]  16. 

which  can  be  solved  for  the  field  inside  V  and  on  S. 

Whereas  the  generation  of  (2)  using  FEM  is  standard,  the  generation  of  (3) 
using  MoM  can  take  many  different  forms.  The  basic  equations  for  generating  (3) 
are  the  electric-field  integral  equation  (EFIE)  given  by 

L(  j)  -  K(M)  =  E*  (5) 

and  the  magnetic-field  integral  equation  (MFIE)  given  by 

K(  j)  -I-  L(M)  =  H*  (6) 

where  J  and  M  are  related  to  the  fields  on  5”  by  J  =  n  x  H  and  M  =  E  x  n, 
respectively,  and  (E‘,H’)  denote  the  incident  fields.  The  operators  L  and  K  axe 
defined  as 

L(X)  =  A  X(r')  +  ^ VV  .  X(r')  G(r,  (7) 

K(X)  =  TY(r)  +  /  X(r')  x  VG(r,  r')dS'  (8) 

V  S 

where  Y  is  related  to  X  by  X  =  n  x  Y  and 

in  which  i?  =  |r-r'|.  The  bar  integral  symbol  is  used  to  denote  the  principal  value 
and  the  parameter  T  is  given  by  T  =  1  — n/47r  where  ft  is  the  solid  angle  subtended 
by  the  observation  point  [19].  For  a  smooth  surface,  ft  =  27r  and  T  =  1/2. 
Equations  (5)  and  (6)  can  be  discretized  by  first  expanding  J  and  M  as 


Ns 

(10) 

t=l 

Ns 

t=l 

(11) 
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where  Ns  denotes  the  total  number  of  edges  on  S  and  gi  denotes  the  RWG  vector 
basis  functions  [20],  which  are  completely  compatible  with  the  vector  basis  func¬ 
tions  for  the  edge  elements.  Substituting  (10)  and  (11)  into  (5)  and  using  gj  as  the 
weighting  function,  we  obtain  the  TE  formulation  (short  for  I  •  E  where  i  denotes 
a  unit  vector  tangential  to  S) 


+  [(2™]{ffs}  =  {6™} 

(12) 

where 

(13) 

Qjf  =  lg>-Uei)dS 

(14) 

6P=  /  g,  E‘(iS. 

J  s 

(15) 

Similarly,  from  (6) 

we  obtain  the  TH  formulation  (short  for  i  •  H) 

1/’™]{£5}  +  =  {t™} 

(16) 

where 

PFi’'  =  l&-MSi)dS  =  Qjf 

(17) 

«5»=/^a-K(g,)d5  =  -pJ^ 

(18) 

6™=  /  gi-ff  (<5. 

*/  s 

(19) 

Alternatively,  we  may  choose  h  x  g,-  as  the  weighting  function  and  obtain  from 
(5)  the  NE  formulation  (short  for  n  x  E) 

(20) 

where 

P<r  =  -  l«‘X»-K{si)dS 

(21) 

/  ”  X  g«  'L(gi) 

J  S 

(22) 

6f^  =  j^nxg,  -E’  dS 

(23) 

and  from  (6)  the  NH  formulation  (short  for  n  x  H) 

(24) 
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where 


/>r=j 

f  n  X  g,  ■  L(gi)  dS  =  Oj® 

s 

(25) 

or=j 

f  nxgi-K(gj)dS  =  -P’^^ 
s 

(26) 

II 

f  nxgi-  H‘  dS. 

(27) 

Equations  (20)  and  (24)  can  also  be  obtained  by  taking  the  cross  product  of  n  with 
(5)  and  (6)  and  then  use  gi  as  the  weighting  function  (That  is  the  reason  we  used 
the  abbreviations  NE  and  NH  for  the  two  equations). 

Theoretically,  any  of  (12),  (16),  (20),  and  (24)  can  be  used  as  (3).  However, 
each  of  them  suffers  from  the  problem  of  interior  resonance  and  fails  to  produce 
accurate  solution  at  and  near  certain  frequencies  corresponding  to  the  resonant 
frequencies  of  the  cavity  formed  by  covering  S  with  a  perfect  electric  or  magnetic 
conductor  and  filling  it  with  the  exterior  medium.  To  eliminate  this  problem,  one 
has  to  combine  an  equation  from  EFIE  to  another  equation  from  MFIE  to  obtain 
a  combined  equation  (that  is,  CFIE)  [21].  For  example,  one  can  combine  (12) 
with  (16)  to  obtain  the  TETH  formulation  or  (12)  with  (24)  to  obtain  the  TENH 
formulation.  One  can  also  combine  (20)  with  (16)  to  obtain  the  NETH  formulation 
or  (20)  with  (24)  to  obtain  the  NENH  formulation  which  is  the  one  employed  in 
[12].  Among  the  four  CFIE  combinations,  TENH  and  NETH  are  used  most  widely. 
However,  it  is  not  clear  which  combination  would  produce  the  most  efficient  and 
accurate  solution. 

Let  us  consider  the  issue  of  efficiency  first.  It  is  known  that  the  FEM  matrices 
in  (4)  are  diagonally  dominant.  Hence,  (4)  would  be  better  conditioned  if  [(5]  is 
diagonally  dominant.  An  analysis  of  the  matrix  property  shows  that  [P^®]  and 
[qJVH]  diagonally  dominant,  [Q^^]  and  [P^^]  are  diagonally  dominant, 

and  [P^^],  [Q^^],  [Q^^],  and  [P^^]  are  least  diagonally  dominant.  These  facts 
can  be  denoted  as 
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From  these,  we  obtain  the  matrix  structure  for  the  TETH  formulation  as 

I  ••• 

[P\Q]  ~  1  I  1  .  (29) 

I 

For  the  TENH  formulation,  we  have 


For  the  NETH  formulation,  we  have 


3  I  0  .  (31) 

■••  I 

Finally,  for  the  NENH  formulation,  we  have 


Considering  the  properties  of  the  FEM  matrices  in  (4),  heuristically,  it  is  ap¬ 
parent  that  the  TENH  formulation  would  produce  the  best  conditioned  matrix  for 
(4),  the  NETH  formulation  would  yield  the  worst  conditioned  matrix,  and  both 
TETH  and  NENH  formulations  have  condition  numbers  between  those  of  TENH 
and  NETH. 

To  verify  the  above  predictions,  we  consider  the  problem  of  plane-wave  scat¬ 
tering  by  a  coated  sphere.  The  coated  sphere  has  a  radius  rj  and  its  conducting 
core  has  a  radius  r^.  The  dielectric  coating  heis  a  relative  permittivity  =  4  and 
a  free-space  permeability  and  its  thickness  is  chosen  large  enough  so  that  there 
is  an  appreciable  tangential  electric  field  on  the  surface.  Equation  (4)  is  solved 
using  the  conjugate  gradient  (CG)  method.  Figure  1  displays  the  residual  norm 
versus  the  number  of  iterations,  from  which  we  see  clearly  that  TENH  converges 
most  quickly,  NETH  has  the  worst  convergence,  and  the  convergence  of  TETH  and 
NENH  lies  between  those  of  TENH  and  NETH.  This  observation  agrees  perfectly 
with  our  earlier  prediction. 
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Next,  we  consider  the  issue  of  accuracy.  Examining  (12)-(27)  carefully,  we  find 
that  in  the  TE  formulation,  where  gj  is  used  as  the  weighting  function,  the  first 
term  in  (8)  has  no  contribution  to  (13)  when  i  =  j,  or  in  other  words,  the  first  term 
in  (8)  is  not  tested.  The  same  observation  can  be  made  for  the  TH  formulation. 
However,  in  the  NE  formulation,  where  h  x  g,-  is  used  as  the  weighting  function, 
the  first  term  in  (7)  cannot  be  tested,  and  thus  has  no  contribution  to  (22)  when 
i  =  j.  The  same  observation  can  be  made  for  the  NH  formulation.  Clearly,  neither 
gj  nor  n  X  gi  forms  a  complete  set  of  weighting  function  for  (5)  or  (6).  Therefore, 
when  gi  or  n  X  g,  is  used  alone,  the  solution  can  become  inaccurate.  Since  all 
the  formulations  described  earlier  (TETH,  TENH,  NETH,  and  NENH)  are  the 
result  of  using  either  g,-  or  n  x  g,  as  the  weighting  function,  their  solutions  can  be 
inaccurate  as  well. 

The  above  analysis  on  accuracy  is  also  verified  by  the  numerical  analysis  of  the 
problem  described  earlier.  Figure  2  shows  the  bistatic  radar  cross  section  (RCS) 
of  the  coated  sphere.  It  is  obvious  that  all  the  four  formulations  have  a  significant 
error  in  their  solutions.  Our  further  numerical  experiments  show  that  such  errors 
cannot  be  reduced  significantly  by  using  finer  discretization.  It  is  interesting  to 
note  that  both  TETH  and  NENH  have  a  similar  error  and  both  TENH  and  NETH 
also  have  a  similar  error.  However,  the  error  in  TETH  and  NENH  is  smaller  than 
that  in  TENH  and  NETH.  We  note  that  this  problem  of  inaccuracy  occurs  only 
when  there  exist  simultaneously  nontrivial  tangential  electric  and  magnetic  fields 
on  the  surface  5;  therefore,  it  disappears  when  one  deals  with  a  bare  conducting 
body  or  a  conducting  body  with  a  very  thin  coating  where  the  tangential  electric 
field  is  very  small.  This  is  why  the  problem  was  not  observed  in  [12].  We  also  note 
that  this  problem  was  not  observed  in  [10],  [11],  and  [14]  because  none  of  them 
employed  the  RWG  functions  as  both  the  expansion  and  weighting  functions. 

To  alleviate  the  inaccuracy  discussed  above,  it  is  clear  that  a  more  complete 
set  of  weighting  functions  has  to  be  used.  A  natural  choice  is  a  combination  of  g,- 
and  n  x  g,.  When  this  is  applied  to  (5),  we  obtain  a  matrix  equation,  which  is 
equivalent  to  the  sum  of  (12)  and  (20)  and  is  referred  to  as  the  TENE  formulation. 
When  this  is  applied  to  (6),  we  obtain  a  matrix  equation,  which  is  equivalent  to 
the  sum  of  (16)  and  (24)  and  is  referred  to  as  the  THNH  formulation.  However, 
since  TENE  comes  from  EFIE  and  THNH  comes  from  MFIE,  both  would  suffer 
from  the  problem  of  interior  resonance.  One  remedy  is  to  combine  TENE  and 
THNH.  A  more  efficient  alternative  is  to  combine  TENE  with  either  the  NH  or 
TH  formulation.  A  simple  analysis  of  matrix  condition  shows  that  among  NH  and 
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TH,  NH  is  a  better  choice  for  the  combination  with  TENE.  Figure  3  shows  the 
result  of  TENENH,  from  which  we  see  that  TENENH  has  a  significantly  better 
accuracy  than  those  in  Fig.  2.  The  remaining  error  in  TENENH  can  be  reduced 
by  using  a  finer  discretization.  The  corresponding  convergence  curves  are  given  in 
Fig.  4. 

The  results  presented  above  are  obtained  at  the  frequency  that  does  not  coin¬ 
cide  with  the  frequency  of  interior  resonance.  To  ensure  the  validity  of  our  analysis, 
we  consider  the  same  coated  sphere  at  a  frequency  of  interior  resonance.  Figure  5 
displays  the  residual  norm  versus  the  number  of  iterations,  from  which  we  observe 
a  similax  convergence  behavior  that  agrees  with  our  prediction.  However,  com¬ 
pared  to  Fig.  1,  the  number  of  iterations  for  TETH  and  NENH  in  this  case  has 
increased  significantly  whereas  that  for  TENH  and  NETH  remains  the  same.  To 
investigate  this  problem  further,  we  recorded  the  number  of  iterations  at  the  fre¬ 
quencies  near  the  frequency  of  interior  resonance  and  the  result  is  given  in  Fig.  6. 
To  our  surprise,  both  TETH  and  NENH  have  a  sharp  peak  at  the  frequency  of  in¬ 
terior  resonance.  This  implies  that  both  TETH  and  NENH  yield  an  ill-conditioned 
matrix  and  still  suffer  from  the  problem  of  interior  resonance,  although  they  are 
derived  from  the  CFIE  formulation.  However,  the  bandwidth  of  the  ill-conditioned 
peaks  is  extremely  narrow  (less  than  1%),  compared  to  those  resulting  from  either 
the  EFIE  or  the  MFIE  (about  10%),  and  this  is  probably  the  reason  that  this 
problem  was  not  detected  before.  The  results  for  the  RCS  are  given  in  Fig.  7. 
As  expected,  both  TETH  and  NENH  yield  a  result  drastically  different  from  the 
exact  solution  whereas  both  TENH  and  NETH  produce  a  stable  result  with  an 
error  similar  to  that  in  Fig.  2.  The  result  obtained  using  TENENH  is  presented  in 
Fig.  8,  from  which  a  good  agreement  is  observed.  The  number  of  iterations  at  the 
frequencies  near  the  frequency  of  interior  resonance  is  also  given  in  Fig.  6,  showing 
a  very  stable  behavior. 

Next,  we  present  several  other  examples  to  demonstrate  the  accuracy  and  ca¬ 
pability  of  the  proposed  formulation  for  other  geometries.  Figure  9  shows  the 
bistatic  RCS  of  a  finite  dielectric  cylinder  and  Fig.  10  displays  the  result  for  a 
dielectric  cube.  All  the  results  are  compared  with  those  obtained  from  MoM  and 
excellent  agreement  is  observed  in  each  case.  We  note  that  the  MoM  solutions 
shown  in  Figs.  9  and  10  are  obtained  from  the  PMCHW  formulation  [22],  which  is 
a  combined-source  integral  equation  (CSIE).  The  PMCHW  formulation  is  known 
to  produce  accurate  solution  [23],  [24];  however,  it  can  be  applied  to  only  homoge¬ 
neous  objects. 
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III.  CONCLUSION 

In  this  paper,  we  studied  in  detail  a  variety  of  formulations  for  the  hybrid  FE-BI 
method  for  calculating  3D  electromagnetic  scattering  by  inhomogeneous  objects. 
We  showed  that  the  efficiency  and  accuracy  of  the  FE-BI  method  depend  highly  on 
the  formulation  and  discretization  of  BIE  used.  We  considered  four  formulations 
(TETH,  TENH,  NETH,  and  NENH)  obtained  from  the  discretization  of  the  CFIE, 
and  we  found  from  analysis  and  verified  numerically  that 

•  TENH  produces  the  best  condition  FE-BI  matrix  equation  and  NETH  produces 

the  worst  conditioned  matrix  equation.  Therefore,  when  an  iterative  solver 
such  as  the  CG  algorithm  is  employed  to  solve  the  matrix  equation,  TENH 
is  the  most  efficient  formulation. 

•  None  of  the  four  formulations  produces  accurate  FE-BI  solution  because  neither 

the  RWG  vector  basis  functions  (g,)  nor  its  cross  product  with  the  unit 
normal  (n  x  gi)  form  a  complete  set  of  weighting  functions  for  EFIE  or 
MFIE  on  a  general  surface  where  nontrivial  equivalent  electric  and  magnetic 
currents  exist  simultaneously. 

•  Both  TETH  and  NENH  suffer  from  the  problem  of  interior  resonance  although 

the  bandwidth  of  the  bad  solution  is  extremely  narrow,  compared  to  those 
resulting  from  EFIE  and  MFIE.  However,  TENH  and  NETH  are  immune  to 
the  problem  of  interior  resonance  although  their  results  are  inaccurate. 

Based  on  the  analysis,  we  proposed  a  formulation  (TENENH)  that  has  a  good 
efficiency  and  a  good  accuracy  and  is  completely  immune  to  the  corruption  of 
interior  resonance.  The  TENE  part  of  this  formulation  is  equivalent  to  testing  the 
pertinent  EFIE  twice,  or  equivelent  to  testing  an  equation  with  N  unknowns  2N 
times,  yielding  2N  equations.  This  is  an  overdetermined  system  whose  solution 
can  be  sought  by  the  least  square  approach.  However,  we  find  that  by  adding  the 
equations,  which  is  equivalent  to  testing  EFIE  by  gi  -1-  n  x  gi,  good  result  is  also 
obtained.  Our  next  step  is  then  to  apply  MLFMA  to  the  proposed  FE-BI  method 
to  enhance  its  capability  to  deal  with  larger  objects. 
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Figure  1:  The  normalized  residual  norm  versus  the  number  of  iterations  in  the  CG 
solution  of  scattering  by  a  coated  sphere. 
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Figure  2:  The  bistatic  RCS  of  a  coated  sphere.  Neither  of  the  four  formulations 
produces  accurate  results. 
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Figure  3:  The  bistatic  RCS  of  a  coated  sphere.  Good  results  are  obtained  using 
the  TENENH  formulation. 


Number  of  Iterations 

Figure  4:  The  normalized  residual  norm  versus  the  number  of  iterations  for  the 
TENENH  formulation. 
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Figure  5:  The  normalized  residual  norm  versus  the  number  of  iterations  in  the  CG 
solution  of  scattering  by  a  coated  sphere  at  the  frequency  of  interior  resonance. 
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Figure  6:  The  number  of  iterations  versus  frequency  (equivalently,  the  size  of  the 
scatterer  in  terms  of  wavelength).  Both  NENH  and  TETH  exhibit  singular  behav¬ 
ior  near  the  frequency  of  interior  resonance  whereas  both  TENH  and  TENENH 
display  a  stable  behavior. 
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Figure  7:  The  bistatic  RCS  of  a  coated  sphere.  Neither  of  the  four  formulations 
produces  accurate  results.  In  particular,  both  NENH  and  TETH  yield  erroneous 
results. 
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Figure  8:  The  bistatic  RCS  of  a  coated  sphere  at  the  frequency  of  interior  reso¬ 
nance.  Again,  good  results  are  obtained  using  the  TENENH  formulation. 
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Figure  9:  The  bistatic  RCS  of  a  finite  dielectric  cylinder  in  the  x-z  plane  for  a 
plane  wave  incident  along  the  2-axis,  (a)  VV-polarization.  (b)  HH-polarization. 
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Figure  10:  The  bistatic  RCS  of  a  dielectric  cube  in  the  x-z  plane  for  a  plane  wave 
incident  along  the  z-axis.  (a)  VV-polarization.  (b)  HH-polarization. 
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Abstract 

A  new  approach  is  proposed  to  reduce  the  reflection  error  of  a  perfectly  matched  layer  (PML)  in 
the  frequency-domain  finite-element  solution  of  electromagnetics  problems.  The  approach  is  based  on 
the  complementary  nature  of  a  PML  backed  by  a  perfect  electric  conductor  (PEC)  and  one  backed  by 
a  perfect  magnetic  conductor  (PMC).  By  averaging  the  solutions  obtained  by  PEC-  and  PMC-backed 
PML,  the  error  is  reduced  by  20log|i?|  dB  where  R  is  the  reflection  coefficient  of  the  PEC-  or  PMC-backed 
PML.  Numerical  results  are  presented  to  demonstrate  the  accuracy  of  the  complementary  PMLs. 


1  Introduction 

The  perfectly  matched  layer  (PML)  has  recently  been  introduced  by  Berenger  [1]  as  a  material  absorbing 
boundary  condition  (ABC)  for  electromagnetics  problems.  An  infinitely  thick  PML  has  no  reflection  for  all 
incident  angles  and  all  frequencies  at  its  interface  with  another  medium.  However,  when  the  PML  is  used 
to  truncate  the  finite  difiference  or  finite  element  solution  domain,  it  must  be  terminated  or  backed  by  a 
perfectly  conducting  or  impedance  surface.  As  a  result,  the  PML  does  not  possess  zero  reflection  for  all 
incident  angles.  Instead,  the  reflection  coefficient  increases  as  the  incident  angle  and  reaches  1  at  grazing. 
To  reduce  the  error  from  this  artificial  reflection,  one  must  either  place  the  PML  at  some  distance  away 
from  the  object  or  increase  the  PML  thickness.  Both  would  result  in  a  large  solution  domain  and  reduce  the 
efficiency  of  the  numerical  solution.  Recently,  we  have  studied  the  optimization  of  the  PML  [2]  and  proposed 
a  combination  of  the  PML  with  ABC  to  reduce  the  artificial  reflection  [3]. 

•The  authors  are  with  the  Center  for  Computational  Electromagnetics,  Department  of  Electrical  and  Computer  Engineering, 
University  of  Illinois  at  Urbanar Champaign,  Urbana,  Illinois  61801-2991  USA. 
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In  this  paper,  we  propose  a  new  approach  to  reduce  the  numerical  error  associated  with  the  artificial 
reflection  of  the  PML.  This  approach  employs  a  pair  of  PMLs  that  have  exactly  the  same  material  composition 
or  parameters,  but  one  is  backed  by  a  perfect  electric  conductor  (PEC)  and  the  other  by  a  perfect  magnetic 
conductor  (PMC).  As  a  result,  the  PEC-backed  PML  has  a  reflection  coefficient  negative  to  that  of  the 
PMC-backed  PML,  or  in  other  words,  the  PEC-  and  PMC-backed  PMLs  are  complementary  to  each  other. 
Assuming  that  the  PEC-backed  PML  has  a  reflection  coefficient  R,  the  numerical  solution  obtained  with  this 
PML  has  an  error  proportional  to  R.  When  the  same  problem  is  solved  using  the  PMC-backed  PML,  the  error 
is  proportional  to  —R.  Averaging  the  two  solutions  eliminates  the  error  proportional  to  R  and  the  remaining 
error  is  proportional  to  R? ,  which  is  due  to  the  double  reflections  of  waves.  The  use  of  the  complementary 
PMLs  is  demonstrated  by  both  two-  and  three-dimensional  waveguide  discontinuity  problems.  We  note 
that  a  similar  idea  was  employed  before  for  the  finite-difference  time-domain  (FDTD)  solution  of  partial 
differential  equations  [4-6]. 

2  Analysis 

To  demonstrate  the  application  of  complementary  PMLs,  consider  a  parallel-plate  waveguide  discontinuity 
problem,  sketched  in  Fig.  1(a).  To  solve  this  problem  using  the  finite-element  method  (FEM)  in  conjunction 
with  the  PML,  we  terminate  both  ends  with  the  PML  and  place  a  magnetic  current  sheet  to  excite  the 
desired  incident  field,  as  illustrated  in  Fig.  1(b).  Using  the  coordinate-stretching  approach  proposed  in  [7] 
and  extended  in  [3],  we  can  show  that  for  a  TM  mode  incidence,  the  magnetic  field  satisfies  the  differential 
equation 
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where  e*  and  are  coordinate  stretching  variables.  For  a  PML  normal  to  the  x  axis,  e*  =  Si(x)  and  =  1, 
and  for  a  PML  normal  to  the  2:  axis,  -  s^(r:)  and  c*  =  1,  where  s*  and  Sz  are  parameters  for  the  PML. 
In  this  case,  there  is  no  PML  normal  to  the  x  axis  and  for  the  PML  normal  to  the  z  axis,  we  choose 


(^)  —  1  ~  i^tnax  (1/ L)^  (2) 

where  Smax  denotes  the  maximum  loss  tangent,  I  is  the  distance  from  the  air/PML  interface,  and  L  is  the 
thickness  of  the  PML. 

To  illustrate  the  accuracy  of  the  complementary  PMLs,  we  first  considered  the  parallel-plate  waveguide 
without  the  discontinuity  and  calculated  the  field  distribution  for  the  TEM  incidence  at  15  GHz.  For  the 
calculation,  L  =  5  mm  and  Smax  =  30//  where  /  denotes  the  frequency  in  GHz.  For  such  a  PML,  the 
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reflection  coeflScient  is  6.3%  (—24  dB),  The  FEM  discretization  is  1  element  per  1  mm.  The  relative  error 
in  the  solution  obtained  with  the  PEC-backed  PML  (PEC-PML)  and  complementary  PMLs  (COM-PML) 
is  given  in  Fig.  2.  The  error  in  the  PMC-PML  solution  is  similar  to  that  in  the  PEC-PML  solution.  As  can 
be  seen,  the  maximum  error  in  the  PEC-PML  solution  is  about  6.6%  and  that  in  the  COM-PML  solution  is 
about  0.33%  (—50  dB),  which  agrees  with  our  earlier  prediction,  that  is,  about  the  square  of  6.6%. 

The  reflection  coefficient  due  to  the  discontinuity  for  the  TEM  incidence  is  given  in  Fig.  3,  where  we  have 
plotted  the  results  obtained  with  PEC-PML,  PMC-PML,  and  COM-PML,  compared  to  the  solution  obtained 
using  the  exact  ABC  [8].  Excellent  agreement  is  observed  between  the  COM-PML  and  exact  ABC  solutions, 
whereas  those  obtained  with  PEC-PML  and  PMC-PML  deviate  from  the  exact  solution  significantly.  We 
note  that  with  the  COM-PML,  we  can  place  the  PML  as  close  as  3  mm  away  from  the  nearest  edge  of  the 
discontinuity  (resulting  in  226  unknowns)  with  an  accuracy  comparable  to  that  obtained  with  a  20-layer 
PEC-backed  PML  placed  10  mm  away  from  the  discontinuity  (resulting  in  710  unknowns)  [2].  Therefore, 
although  the  COM-PML  requires  solving  the  problem  twice,  it  is  still  more  efficient  than  the  PEC-PML 
because  of  the  reduced  number  of  unknowns. 

Next,  we  consider  a  three-dimensional  waveguide  discontinuity  problem,  illustrated  in  Fig.  4.  To  numer¬ 
ically  solve  this  problem,  we  again  terminate  both  ends  of  the  waveguide  with  the  PML  and  introduce  an 
electric  current  to  excite  the  desired  incident  field.  The  differential  equation  for  the  electric  field  is  given  by 

[3] 


VeX 


—  (Ve  X  E) 
Pr 


-  klcrE  =  -jkoZoJ 


(3) 


where 


^  .  1  5  1  .  1  5  1  .  1  a  1 

Ve  =  — 7=  +  — ;=  +  — =  (4) 

y/^OXyJei  y/^OVy/^  y/^dZy/tz 

in  which  Ca?,  Cy,  and  are  defined  in  the  same  way  as  above. 

Again,  to  illustrate  the  accuracy  of  the  COM-PML,  we  calculated  the  field  in  the  empty  waveguide  due  to 
the  TEio  incidence  at  10  GHz.  For  the  calculation,  L  =  6  mm  and  S^ax  =  30//,  and  the  FEM  discretization 
is  1  element  per  2  mm.  The  relative  error  is  given  in  Fig.  5,  where  the  maximum  error  in  the  PEC-PML 
solution  is  about  15%,  in  contrast  to  1.4%  in  the  COM-PML  solution. 

The  reflection  coefficient  due  to  the  discontinuity  for  the  TEio  incidence  is  shown  in  Fig.  6.  For  the 
calculation,  the  PML  is  placed  only  2  elements  away  from  the  nearest  edge  of  the  dielectric  obstacle.  The 
results  are  compared  to  the  data  obtained  using  the  method  of  orthogonal  expansions  [9].  Excellent  agreement 
is  observed  between  the  COM-PML  and  reference  data. 
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3  Conclusion 


A  new  and  simple  approach  was  proposed  to  significantly  reduce  the  PML  reflection  error  in  the  finite 
element  solution  of  electromagnetics  problems.  The  approach  was  validated  by  numerical  examples  for  both 
two-  and  three-dimensional  waveguide  discontinuity  problems.  It  was  shown  that  the  reflection  error  can  be 
reduced  from  201og|J?|  dB  to  401og|i?|  dB. 
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FIGURE  CAPTIONS 


Figure  1.  (a)  Parallel-plate  waveguide  discontinuity  problem,  (b)  PML  truncation. 

Figure  2.  Relative  error  in  the  field  of  an  empty  parallel-plate  waveguide  when  the  waveguide  is  terminated 
by  (a)  PEC-PML  and  (b)  COM-PML. 

Figure  3.  Reflection  coefficient  for  the  dominant-mode  incidence. 

Figure  4.  Rectangular  waveguide  loaded  with  a  dielectric  obstacle  having  €r  =  6.  For  numerical  simulation, 
the  waveguide  is  terminated  with  a  PML. 

Figure  5.  Relative  error  in  the  field  of  an  empty  rectangular  waveguide  when  the  waveguide  is  terminated 
by  (a)  PEC-PML  and  (b)  COM-PML. 

Figure  6.  Reflection  coefficient  for  different  dielectric  lengths,  (a)  d  =  1.2  cm.  (b)  d  =  0.8  cm. 
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Figure  6.  Reflection  coefficient  for  different  dielectric  lengths,  (a)  d  =  1.2  cm.  (b) 
d  =  0.8  cm. 
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Simple  and  Efficient  Computation  of  Electromagnetic  Fields 
in  Arbitrarily-Shaped,  Inhomogeneous  Dielectric 
Bodies  Using  Transpose-Free  QMR  and  FFT 

C,  F.  Wang  and  J.  M.  Jin 

Center  for  Computational  Electromagnetics 
Department  of  Electrical  and  Computer  Engineering 
University  of  Illinois  at  Urbana-Champaign 
Urbana,  Illinois  61801-2991 

Abstract 

A  simple  and  efficient  numerical  method  is  presented  for  computing  electromagnetic 
fields  in  three-dimensional  inhomogeneous  dielectric  bodies.  The  method  employs  a  two- 
stage  discretization  to  convert  an  integro-differential  equation  into  an  implicit  system  of 
linear  algebraic  equations.  This  discrete  system  is  then  solved  using  a  transpose-free  quasi- 
minimal  residual  (TFQMR)  algorithm,  which  avoids  the  complicated  calculation  of  the 
multiplication  between  the  transpose  of  the  system  matrix  and  a  vector.  The  simple  mul¬ 
tiplication  between  the  system  matrix  and  a  vector  required  in  the  TFQMR  algorithm  is 
calculated  eflSciently  using  only  six  fast  Fourier  transforms  (FFT).  Numerical  results  for 
strongly  inhomogenous  and  lossy  spheres  show  that  the  method  has  a  stable  convergence 
behavior  and  excellent  numerical  performance. 

I.  Introduction 

Efficient  computation  of  electromagnetic  fields  in  arbitrarily-shaped,  inhomogeneous 
dielectric  bodies  in  a  three-dimensional  (3D)  space  plays  an  important  role  in  many  appli¬ 
cations  such  as  nondestructive  testing,  microwave  imaging,  scattering  control,  target  identi¬ 
fication,  electromagnetic  hyperthermia,  and  magnetic  resonance  imaging.  The  well-known 
method  of  moments  (MoM)  [l]-[3]  is  one  of  the  popular  methods  for  this  computation. 
In  this  method,  an  integro-differential  equation  is  first  formulated  in  terms  of  volumetric 
equivalent  current  that  accounts  for  the  effect  of  the  permittivity  and  conductivity  of  an 
inhomogeneous  body.  This  integro-differential  equation  is  then  discretized  using  mostly 
Galerkin’s  procedure.  The  discretization  results  in  a  matrix  equation  with  a  very  large 
number  of  unknowns,  whose  solution  using  a  direct  solver,  such  as  Gaussian  elimination 
and  LU  decomposition  method,  is  basically  impractical,  because  a  direct  solver  has  a  mem¬ 
ory  requirement  of  0{N^)  and  computational  complexity  of  0{N^),  where  N  denotes  the 


92 


number  of  unknowns.  This  difficulty  can  be  circumvented  by  solving  the  matrix  equation 
using  an  iterative  solver  and  in  each  iteration  the  required  matrix-vector  multiplication  is 
evaluated  using  the  fast  Fourier  transform  (FFT)  [4].  In  the  past,  the  conjugate  gradient 
(CG)  and  the  biconjugate  gradient  (BCG)  methods  have  been  employed  as  such  an  itera¬ 
tive  solver,  and  the  resultant  methods  are  often  referred  to  as  the  CG-FFT  and  BCG-FFT 
methods  [5]-[7].  The  use  of  the  so-called  CG-FFT  or  BCG-FFT  method  reduces  the  mem¬ 
ory  requirement  to  0{N)  and  computational  complexity  to  0{NiterN  log  N),  where  Nuer 
denotes  the  number  of  CG  or  BCG  iterations. 

There  is  a  large  body  of  literature  on  the  CG-FFT  and  BCG-FFT  methods  for  a 
variety  of  electromagnetics  problems  and  it  is  not  our  intention  to  review  it  here.  Instead, 
we  shall  focus  on  those  for  3D  volumetric  material  problems.  The  first  application  of 
the  CG-FFT  method  to  such  problems  can  be  found  in  the  analysis  of  the  absorption  of 
electromagnetic  power  by  human  bodies  [5].  However,  the  use  of  pulse  basis  functions 
yielded  slow  convergence  and  poor  results  when  dealing  with  materials  with  high  dielectric 
contrast.  Better  formulations  were  later  proposed  [8]-[12],  and  most  used  mixed-order 
(linear  in  one  direction  and  constant  in  the  other  two  directions)  basis  functions.  Among 
these,  the  methods  proposed  by  Zwamborn  and  van  den  Berg  [10]  and  Gan  and  Chew  [11] 
are  the  most  accurate  for  materials  with  high  dielectric  contrast.  In  both  methods,  one  is 
required  to  calculate  within  each  iteration  the  multiplication  between  the  transpose  of  the 
system  matrix  and  a  vector,  in  addition  to  that  between  the  system  matrix  and  a  vector, 
resulting  in  at  least  12  FFTs  per  iteration.  Furthermore,  the  multiplication  between  the 
transpose  of  the  system  matrix  and  a  vector  is  found  to  be  more  complicated  than  that 
between  the  system  matrix  and  a  vector. 

In  this  paper,  we  present  an  alternative  and  more  efficient  method  for  computing  elec¬ 
tromagnetic  fields  in  arbitrarily-shaped,  inhomogeneous  dielectric  bodies.  In  this  method, 
a  transpose-free  quasi-minimal  residual  (TFQMR)  algorithm  [13]  is  employed  to  avoid  the 
complicated  multiplication  between  the  transpose  of  the  system  matrix  and  a  vector,  result¬ 
ing  in  a  much  simpler  computer  implementation.  Moreover,  the  number  of  FFTs  is  reduced 
to  only  six  per  iteration.  It  is  observed  that  the  TFQMR-FFT  method  yields  excellent 
results  even  for  highly  inhomogeneous  dielectric  objects. 

II.  Formulation 

Consider  the  problem  of  scattering  by  a  lossy  inhomogeneous  dielectric  object  with  a 
complex  permittivity 

/X  /X  .o’(x)  ,  . 

e(x)  =  €,.(x)€o-j-^  (1) 

u 

where  Cr  denotes  the  relative  permittivity  and  a  denotes  the  electric  conductivity  of  the 
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object,  which  is  in  a  free-space  having  a  permittivity  cq-  The  incident  electric  field  is 
denoted  as  .  The  scattering  problem  can  be  formulated  as  the 

following  domain  integral  equation  over  the  object  domain  V ; 

E‘”(x)  =  ^-(l=S  +  VV-)A(x),  xeV  (2) 

where  fco  =  wy'coMo  and 

A(x)  =  -  /"  G(x  -  x')x(x')D(x')dx'  (3) 

Co  Jw 

with 


v(x)  = 

^  e(x)  ’ 


G(x  -  x')  = 


exp(-j/:o|x-  x'l) 


47r|x  —  x'l 

To  discretize  this  equation,  we  place  the  object  in  a  uniform  mesh  with  grid  widths  of 
Axi,  Ax2)  and  Ax3  in  the  xi,  X2,  and  X3  directions,  respectively.  Therefore,  the  object  is 
modelled  approximately  as  a  collection  of  small  grids.  The  center  of  each  grid  is  denoted 
as  'x.m,n,P  =  {{M  —  5)Axi,  {N  —  5)Ax2,  {P  -  5)Ax3}  and  within  each  grid  the  complex 
permittivity  is  assumed  to  be  constant  with  value  €m,n,p  —  ^{^M,N,p)' 

To  convert  (2)  into  a  matrix  equation,  we  expand  the  generalized  electric  flux  density 
and  the  electric-contrast  vector  potential  as 

D(x)  =  XI  X  €  V  (4) 


q=lI,J,K 

A(x)  =  XI  H  X  €  V 

q=lI,J,K 


(5) 


where  ^/^j,/f(x),  '5'/^j,k'(x),  are  vector  volumetric  rooftop  functions  in  Xi,  X2 

and  X3  directions,  respectively  [10],  [11].  We  then  apply  the  Galerkin’s  testing  formulation 
to  (2)  and  obtain 

('*'mV,p(x))E’"'^(x))  =  -  kQ{‘9^^f^p{x),A{x)) 

+<V-^&UpW*'^-A(x)>  (6) 

for  g  =  1, 2, 3,  where  (•)  denotes  the  inner  product  of  two  vector  functions.  Substituting  (4) 
and  (5)  into  (6),  we  obtain  the  following  weak  form  of  domain  integral  equation 


P;I,J,K  ~  ‘^M,N,P-,I,J,K. 


(p.s) 


]{AUk\ 


(7) 


where 


tnc,(p)  _  ,yjy(p)  pincv  (p,?)  _  /^(p)  \ 

^M,N,P  -  \^M,N,P''  ^  /’  ^M,N,PiI,J,K  “  \^M,N,P^ 

„(P.9)  _  /m(p)  q/(9)  \  ,,,(P.9)  _  /yj  .  iT,(p)  V  .  \ 
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The  relationship  between  d!j^,N,P  ^M,N,P  found  by  substituting  (4)  and  (5)  into 

(3),  yielding 

^M,N,P=^^  S  Gm-M',N-N',P-P'Xm\N',P'^^M',N',P' 

M',N>,P' 

=  AVDFT-'  {DFT{Gm.n.f)  ■  />F3’{xg„,p<4!„,p})  (») 

where  AV  =  AiiAxaAss.  Substituting  (8)  into  (7),  we  obtain  a  system  of  linear  algebraic 
equations,  which  can  be  symbolically  written  as 

e’""  =  Ld.  (9) 

The  formulation  described  above  was  first  proposed  by  Zwamborn  and  van  den  Berg 
[10].  Its  major  advantage  is  the  simplicity  in  treating  the  singularity  of  the  integrals  in 
(2)  and,  more  important,  in  calculating  the  right-hand  side  of  (7),  which  is  accomplished 
through  two  stages.  The  first  stage  is  to  calculate  ^m,n,p  (®)  second  stage 

is  to  substitute  it  into  (7).  Note  that  the  matrices  implied  in  (7)  are  sparse  matrices  and 
their  product  with  a  vector  can  be  evaluated  with  0{Nt)  operations,  where  Nj  denotes  the 
number  of  unknowns.  Although  the  matrix  implied  in  (8)  is  a  dense  matrix,  the  computation 
of  its  product  with  a  vector  can  be  evaluated  with  0{N’i'\ogN']:)  operations  with  the  aid  of 
the  FFT. 


III.  TFQMR-FFT  Iterative  Algorithm 

Once  (9)  is  formulated,  its  solution  yields  a  numerical  solution  to  the  original  problem. 
However,  since  the  number  of  unknowns  in  (9)  is  usually  very  large,  its  solution  using 
a  direct  solver,  such  as  Gaussian  elimination  and  LU  decomposition  method,  is  basically 
impractical,  because  a  direct  solver  has  a  memory  requirement  of  0{Nj)  and  computational 
complexity  of  0{N^).  This  difficulty  can  be  circumvented  by  solving  (9)  using  an  iterative 
solver  and  in  each  iteration  the  required  matrix-by-vector  product  is  evaluated  using  the 
FFT,  as  pointed  out  earlier.  In  the  past,  the  CG  and  BCG  methods  have  been  employed 
as  such  an  iterative  solver,  and  the  resultant  methods  are  referred  to  as  the  CG-FFT  and 
BCG-FFT  methods  [5]-[12].  The  use  of  these  methods  reduces  the  memory  requirement 
to  0{Nt)  and  computational  complexity  to  ©(AjterArlogiVT),  where  Niter  denotes  the 
number  of  CG  or  BCG  iterations.  However,  both  CG-FFT  and  BCG-FFT  algorithms 
require  the  calculation  of  a  matrix-by-vector  product  with  the  conjugate  transpose  of  the 
system  matrix,  which  is  not  an  easy  task  since  the  system  matrix  in  (9)  is  nonsymmetric. 
Furthermore,  the  CG  method  has  a  problem  of  slow  convergence  although  it  converges 
monotonically,  and  the  BCG  method  does  not  guarantee  convergence  although  it  usually 


95 


converges  quickly.  Although  this  problem  can  be  alleviated  by  using  the  quasi-minimal 
residual  (QMR)  method  [14],  which  converges  monotonically  with  a  convergence  rate  similar 
to  the  BCG  method,  the  QMR  still  requires  the  calculation  of  a  matrix-by-vector  product 
with  the  conjugate  transpose  of  the  system  matrix.  Here,  we  consider  other  alternatives. 

There  are  four  algorithms  that  do  not  require  the  calculation  of  the  transpose  of  the 
system  matrix.  The  first  one  is  the  conjugate  gradient  squared  (CGS)  algorithm  [15],  which 
is  the  transpose-free  variant  of  the  BCG  algorithm.  However,  like  the  BCG  method,  it  also 
exhibits  a  rather  irregular  convergence  behavior  with  wild  oscillations  in  residual  norm  and 
does  not  guarantee  convergence.  The  second  method  is  the  BCG  stabilized  (BCGSTAB) 
algorithm  [16],  which  uses  local  steepest  descent  steps  to  obtain  a  more  smoothly  convergent 
CGS-like  process.  While  this  algorithm  seems  to  work  well  in  many  cases,  it  still  exhibits 
the  irregular  convergence  behavior  for  some  difficult  problems.  Also,  its  convergence  is 
considerably  slower  than  the  CGS  algorithm.  The  third  method  is  the  transpose-free  QMR 
(TFQMR)  method  [13].  This  algorithm  can  be  implemented  easily  by  changing  only  a  few 
lines  in  the  standard  CGS  algorithm.  However,  unlike  the  CGS  algorithm,  the  iterations 
of  the  TFQMR  algorithm  are  characterized  by  a  quasi-minimization  of  the  residual  norm. 
This  leads  to  smooth  convergence  with  a  convergence  rate  similar  to  the  CGS  algorithm. 
The  TFQMR  algorithm  can  be  considered  as  a  new  version  of  the  CGS  algorithm  which 
“quasi-minimizes”  the  residual  in  the  space  spanned  by  the  vectors  generated  by  the  CGS 
iterations.  Recently,  a  QMR  variant  of  the  BCGSTAB  algorithm  (QMRBCGSTAB)  [17] 
is  proposed.  Our  experimental  calculation  shows,  however,  that  its  convergence  can  be 
slower  than  the  TFQMR  for  our  problems.  After  a  comprehensive  comparison,  the  TFQMR 
algorithm  is  chosen  for  this  work. 

The  TFQMR  algorithm  for  solving  (9)  is  given  as  follows  [13],  [18]: 

1.  Compute  Wo  =  uo  =  To  =  e*"'^  -  Ldo,  vq  =  Luo,  go  =  0; 

2.  ro  =  ||ro||,0o  =  »?o  =  0 

3.  Choose  Tq  such  that  po  =  (rj,  to)  ^  0. 

4.  For  m  =  0, 1, 2,.. . ,  until  convergence  Do: 

5.  If  m  is  even  then 

6.  0!m+l  —  ®»n  —  pm/ 

7 .  Uto+1  =  Um  —  Oim'Vm 

8.  Endif 

9.  w,„+i  =  w,n  - 

10.  gm-l-l  —  'Irn  "b  {.^m/^m}Vm^m 

11.  Om+l  =  l|w„+ill2/r„i;Cm+l  =  (l  +  ^m+l)~* 

12.  Tm+l  —  ^m+l  ~ 

13.  —  dm  "b  *l»n+l8m+l 
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14.  If  m  is  odd  then 

15.  Pm+l  =  (WTO+l)ro);/?m-l  =  Pm+l/Pm-\ 

16.  W^.|_j  -f"  Pm—\^7n 

17.  =  Lu^+l  + 

18.  Endif 

19.  EndDo 

The  residual  norm  of  the  approximate  solution  is  given  by  j|rm||  <  (m+  1)2  r^. 
From  the  algorithm,  it  is  easy  to  see  that  each  odd  iteration  requires  two  matrix-by-vector 
products  and  each  even  iteration  does  not  require  a  matrix-by-vector  product.  Therefore, 
on  average  the  TFQMR  algorithm  requires  only  one  matrix-by-vector  product,  which  can 
be  calculated  using  six  FFTs. 

IV.  Numerical  Results 

To  demonstrate  the  accuracy  of  the  TFQMR-FFT  algorithm,  we  analyze  the  scattering 
of  a  plane  wave  from  two  layered  dielectric  spheres  and  compare  the  results  with  the  Mie 
series  solution.  In  all  simulations,  we  assume  that  the  incident  plane  wave  is  polarized  in  the 
Xi  direction  and  propagates  in  the  X3  direction.  The  amplitude  of  the  incident  electric  field 
is  1  V/m.  The  first  sphere  has  two  layers,  whose  inner  layer  has  a  radius  oi  =  0.075  m  and 
€ir  =  72.0  -jT61. 779  and  the  outer  layer  has  a  radius  02  =  0.15  m  and  e2r  =  7.5  -y8.9877. 
The  frequency  is  100  MHz.  The  second  sphere  has  three  layers,  whose  inner  layer  has  a 
radius  ai  =  0.3333Ao  and  6ir  =  1.2,  the  middle  layer  has  a  radius  =  0.6667Ao  and 
€2r  =  2.0,  and  the  outer  layer  has  a  radius  03  =  Ao  and  €3^  =  2.4.  Figure  1  shows 
the  field  in  the  two-layer  dielectric  sphere  and  Fig.  2  gives  the  results  for  the  three-layer 
sphere.  The  convergence  criterion  for  these  results  is  rss  =  ||r,„||/||ro||  <  10“^.  Excellent 
agreement  is  observed  between  the  exact  solution  and  the  numerical  results  obtained  using 
31  X  31  X  31  grids.  The  relative  error  vs.  the  number  of  iterations  is  given  in  Fig.  3  for  both 
cases  for  different  grid  sizes.  As  can  be  seen,  the  TFQMR-FFT  algorithm  exhibits  a  stable 
convergence  behavior.  The  computation  time  needed  to  evaluate  one  iteration,  the  total 
number  of  unknowns,  and  the  required  computer  storage  are  given  in  Table  1. 

To  demonstrate  the  efficiency  of  the  algorithm,  we  consider  again  the  problem  illustrated 
in  Fig.  1.  As  can  be  seen  in  Fig.  3(a),  using  the  TFQMR-FFT  algorithm  with  31  x  31  x  31 
grids,  it  takes  112  iterations  to  reduce  rss  below  10~^.  Since  each  iteration  requires  six 
FFTs,  the  total  number  of  FFTs  is  672.  This  problem  was  also  treated  using  the  CG-FFT 
algorithm  in  [10]  and  the  BCG-FFT  algorithm  in  [11].  For  the  same  grid  size  and  accuracy, 
the  CG-FFT  algorithm  takes  about  360  iterations  and,  since  each  iteration  requires  12 
FFTs,  the  total  number  of  FFTs  is  about  4320,  which  is  6.4  times  that  of  the  TFQMR- 
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Table  1:  Computation  Time  and  Storage  on  DEC  Alpha 


Mesh  size 

FFT  size 

Number  of 

unknowns 

CPU-time 

per  iteration 

Computer 

storage 

15  X  15  X  15 

32  X  32  X  32 

10800 

0.74  sec 

2.2  Mb 

31  X  31  X  31 

64  X  64  X  64 

92256 

9.5  sec 

16  Mb 

63  X  63  X  63 

128  X  128  X  128 

762048 

152  sec 

105  Mb 

FFT  algorithm.  When  the  BCG-FFT  algorithm  is  used,  it  takes  only  54  iterations  and,  since 
each  iteration  requires  18  FFTs,  the  total  number  of  FFTs  is  972,  which  is  1.4  times  that 
of  the  TFQMR-FFT  algorithm.  Furthermore,  the  BCG-FFT  algorithm  has  an  irregular 
convergence  behavior  and  does  not  guarantee  convergence. 

Finally,  to  demonstrate  the  capability  of  the  TFQMR-FFT  algorithm  to  treat  strongly 
inhomogeneous  dielectric  object,  we  consider  the  plane  wave  scattering  by  a  human  head. 
The  construction  of  the  electromagnetic  model  of  the  head  is  discussed  in  [19]  and  the 
material  property  of  the  tissues  of  the  head  is  given  in  [20].  The  incident  wave  propagates 
in  the  —xs  direction  (from  top)  and  the  incident  electric  field  is  polarized  in  the  Si  direction 
(from  the  left  ear  to  the  right  ear).  The  incident  electric  field  has  an  amplitude  of  1  V/m 
and  the  frequencies  considered  are  64  MHz  and  256  MHz.  The  results  are  presented  in  the 
form  of  spatial  absorption  rate  (SAR)  defined  as  SAR  =  a\W[^ /2p,  where  p  denotes  the 
density.  Figures  4  and  5  show  the  SAR  in  the  axial,  sagittal,  and  coronal  planes  at  the 
two  frequencies.  Figure  6  shows  the  relative  error  vs.  the  number  of  iterations  for  the  two 
frequencies.  Again,  the  TFQMR-FFT  algorithm  exhibits  a  stable  convergence  behavior. 


VI.  Conclusion 

This  paper  presented  a  TFQMR-FFT  algorithm  for  computing  electromagnetic  fields  in 
a  3D  arbitrarily-shaped,  inhomogeneous  dielectric  body.  It  is  observed  that  this  algorithm 
yields  excellent  results  and  exhibits  a  very  stable  convergence  behavior.  Because  of  the 
use  of  the  TFQMR  method,  the  algorithm  avoids  the  computation  of  the  multiplication 
between  the  transpose  of  the  system  matrix  and  a  vector,  which  is  required  in  both  the 
CG-FFT  and  BCG-FFT  methods.  As  a  result,  the  programming  complexity  is  greatly 
reduced.  Furthermore,  since  on  average  the  TFQMR  method  requires  only  one  matrix-by¬ 
vector  multiplication,  which  can  be  evaluated  using  six  FFTs,  the  TFQMR-FFT  algorithm 
is  more  efficient  than  the  currently  available  CG-FFT  and  BCG-FFT  methods. 
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Figure  Captions 


Fig.  1  The  magnitude  of  the  total  electric  field  inside  a  two-layer  dielectric  sphere  along 
the  xi,  X2,  and  X3  axes.  The  inner  layer  has  a  radius  oi  =  0.075  m  and  eir  = 
72.0  -  y  161.779  and  the  outer  layer  has  a  radius  02  =  0.15  m  and  62r  =  7.5  -  y8.9877 
and  the  frequency  is  100  MHz.  The  solid  line  is  from  the  Mie  series  solution  and  the 
dash-dot  line  is  from  the  TFQMR-FFT  solution  with  grids  31  x  31  x  31. 

Fig.  2  The  magnitude  of  the  total  electric  field  inside  a  three-layer  dielectric  sphere  along 
the  xi,  X2,  and  X3  axes.  The  inner  layer  has  a  radius  oi  =  0.3333Ao  and  eu  =  1.2, 
the, middle  layer  has  a  radius  02  =  0.6667Ao  and  e2r  =  2.0,  and  the  outer  layer  has  a 
radius  az  =  Ao  and  €3^  =  2.4.  The  solid  line  is  from  the  Mie  series  solution  and  the 
dash-dot  line  is  from  the  TFQMR-FFT  solution  with  grids  31  x  31  x  31. 

Fig.  3  The  relative  error  vs.  the  number  of  iterations  for  different  grid  sizes,  (a)  For  the 
case  of  the  two-layer  sphere,  (b)  For  the  case  of  the  three-layer  sphere. 

Fig.  4  SAR  (W/kg)  in  the  axial,  sagittal,  and  coronal  slices  at  64  MHz  for  a  uniform  plane 
wave  excitation  polarized  in  the  xj  direction  and  propagating  in  the  — X3  direction 
using  63  X  63  X  63  grids. 

Fig.  5  SAR  (W/kg)  in  the  axial,  sagittal,  and  coronal  slices  at  256  MHz  for  a  uniform  plane 
wave  excitation  polarized  in  the  xi  direction  and  propagating  in  the  —  X3  direction 
using  63  X  63  X  63  grids. 

Fig.  6  The  relative  error  vs.  the  number  of  iterations  for  different  frequencies  with  63  x 
63  X  63  grids,  (a)  For  the  case  of  64  MHz.  (b)  For  the  case  of  256  MHz. 


Figure  1:  The  magnitude  of  the  total  electric  field  inside  a  two-layer  dielectric  sphere  along 
the  *1,  X2)  and  xz  axes.  The  inner  layer  has  a  radius  ci  =  0.075  m  and  =  72.0-jT61.779 
and  the  outer  layer  has  a  radius  02  =  0.15  m  and  62r  =  7.5  —  j8.9877  and  the  frequency  is 
100  MHz.  The  solid  line  is  from  the  Mie  series  solution  and  the  dash-dot  line  is  from  the 
TFQMR-FFT  solution  with  grids  31  x  31  x  31. 


Figure  2:  The  magnitude  of  the  total  electric  field  inside  a  three-layer  dielectric  sphere 
along  the  Xi,  X2,  and  X3  axes.  The  inner  layer  has  a  radius  ci  =  0.3333 Aq  and  eir  =  1.2, 
the  middle  layer  has  a  radius  02  =  0.6667Ao  and  €2r  =  2.0,  and  the  outer  layer  has  a  radius 
as  =  Aq  and  (sr  =  2.4.  The  solid  line  is  from  the  Mie  series  solution  and  the  dash-dot  line 
is  from  the  TFQMR-FFT  solution  with  grids  31x31x31. 


102 


Figure  3:  The  relative  error  vs.  the  number  of  iterations  for  different  grid  sizes,  (a)  For  the 
case  of  the  two-layer  sphere,  (b)  For  the  case  of  the  three-layer  sphere. 


Figure  4:  SAR  (W/kg)  in  the  axial,  sagittal,  and  coronal  slices  at  64  MHz  for  a  uniform 
plane  wave  excitation  polarized  in  the  xi  direction  and  propagating  in  the  —x^  direction 
using  63  X  63  X  63  grids. 
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Figure  5:  SAR  (W/kg)  in  the  axial,  sagittal,  and  coronal  slices  at  256  MHz  for  a  uniform 
plane  wave  excitation  polarized  in  the  xi  direction  and  propagating  in  the  —0:3  direction 
using  63  X  63  X  63  grids. 


Figure  6:  The  relative  error  vs.  the  number  of  iterations  for  different  frequencies  with 
63  X  63  X  63  grids,  (a)  For  the  case  of  64  MHz.  (b)  For  the  case  of  256  MHz. 


104 


